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A computer program called Transonic Navier Stokes (TNS) has been devel- 
oped which solves the Euler/Navier-Stokes equations around wings using a zonal 
grid approach. In the present zonal scheme, the physical domain of interest is di- 
vided into several subdomains called “zones” and the governing equations are solved 
interactively. The advantages of the Zonal Grid approach are as follows: l) the grid 
for any subdomain can be generated easily, 2) grids can be, in a sense, adapted to 
the solution, 3) different equation sets can be used in different zones, and 4) this 
approach allows for a convenient data base organization scheme. Using this code, 
separated flows on a NACA 0012 section wing and on the NASA Ames WING 
C have been computed. First, the effects of turbulence and artificial dissipation 
models incorporated into the code are assessed by comparing the TNS results with 
other CFD codes and experiments. Then a series of flow cases are described where 
data are available. The computed results, including cases with shock-induced sep- 
aration, are in good agreement with experimental data. Finally, some “futuristic” 
cases arc presented to demonstrate the abilities of the code for massively separated 
cases which do not have experimental data. 
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CHAPTER 1 


INTRODUCTION 


1.1 Background 

Ever since the introduction of the Navier-Stokes equations in the early 1800’s 
by Navier [l] and Stokes [2], we have been confronted with a very challenging 
situation. Because of the nonlinear and complex nature of these equations, the 
number of closed form analytical solutions has been extremely limited. Hence, 
the Navier-Stokes equations were largely unused for several decades, and did not 
contribute to the early development of airplanes. 

The first major breakthrough in the solution of these equations was intro- 
duced by Ludwig Prandtl [3] in 1904. Prandtl brilliantly reasoned that for high 
Reynolds number flows, viscous effects would be confined to a thin layer along a 
solid surface, called the “Boundary Layer.” We owe much of our understanding 
of the viscous flow phenomenon to this concept. Although this region occupies 
only a small portion of a total flow field, all momentum, heat and mass transfer 
to or from the surface must take place through this layer. Outside the bound- 
ary layer, fluid behaves like an inviscid flow. This leads us to a very important 
conclusion that separate equation sets can be solved in separate domains. There- 
fore, aeronautical problems can be approached from different avenues; i.e., using 
different equation sets, provided that limitations are kept in proper perspective. 
The aforementioned “different equation sets” are subsets of the Navier-Stokes 
equations. If the viscous terms are neglected, the resulting equations describe 
the motion of an inviscid but nonlinear, compressible and rotational fluid: the 
so-called “Euler equations.” If the rotationality is ignored, we get nonlinear, com- 
pressible and irrotational flow; i.e., the full potential equation which can support 
only weak shocks. Further simplification leads to Laplace’s equation which de- 
scribes an incompressible, inviscid, irrotational and steady flow. 

Introduction of digital computers brought a great change to the aeronau- 
tical sciences. They opened new avenues and enabled researchers to tackle the 
problem of solving the above equation sets. With the advent of computers, a 
new discipline of aerodynamics was born, i.e., Computational Fluid Dynamics 
(CFD), and researchers attacked problems which were previously unsolvable. 
But it was soon realized that compromises had to be made because computer 
time and memory limit the range of physical and geometrical scales that can be 
resolved. The more complex the geometry, the less physics we can extract. In 
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terms of numerical techniques, the choice between a complicated set of equations 
with simple geometry and a sophisticated geometry with simpler equations had 
to be made. 

Basically there are three compelling motivations for vigorously developing 
Computational Aerodynamics (Chapman |4]). The first one is to provide new 
technological capabilities which the wind tunnels cannot offer. For example, 
the Reynolds number associated with aircraft flight, the flowfield temperatures 
around atmospheric entry vehicles like the space shuttle, hot jet plumes, or store 
ejections cannot be easily simulated by wind-tunnels. A second compelling mo- 
tivation concerns energy conservation. Wind tunnels consume large amounts of 
energy with increasing unit costs every year, whereas computers do not consume 
as much energy. The third motivation can be explained in terms of economics. 
Chapman [4] indicated that the net cost to conduct a given numerical simula- 
tion with a fixed algorithm has decreased ten times every eight years because 
computer speed has increased at a much greater rate than computer cost. 

Developments of computational aerodynamics have been paced by computer 
memory and speed. Naturally, work started with the most simplified equation 
sets and configurations and has progressed to more sophisticated geometries with 
better approximation of the flow physics. Chapman [4] presents the historical 
progress of Computational Aerodynamics in four stages, each of which is a suc- 
cessively more refined approximation to the full Navier-Stokes equations. Stage I 
is the linearized inviscid approximation or Laplace’s equation. Computational 
methods which solve this governing equation are called “Panel Methods” and 
have been very useful in the design process. The linearized inviscid approxima- 
tion contains only 3 terms; whereas the full Navier-Stokes equations represent- 
ing conservation of mass, momentum and energy contain altogether 60 partial- 
derivative terms when written out in three Cartesian coordinates. Panel methods 
provide realistic determinations of pressure distribution, of lift and side forces 
in attached subsonic and supersonic flow using small disturbance theory. For an 
extensive list of references and survey articles the reader is referred to [4], [5]. 

The second stage of approximation is called “nonlinear inviscid” and here 
only the viscous terms are neglected, but in its most complex form, 27 out of 
60 partial derivative terms are retained from the full Navier-Stokes equations. 
The transonic small-disturbance equation, the full potential equation and the 
Euler equations fall into this category. Nonlinear inviscid computations for tran- 
sonic flow are extensive design tools for the aircraft industry today. For more 
information please refer to [4]. 

The third and fourth stage of Computational Aerodynamics involve the 
Reynolds-averaged Navier-Stokes equations and the Large Eddy Simulation(LES) 
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respectively. The Reynolds-averaged Navier-Stokes equations are obtained by re- 
placing the dependent variables in the original Navier-Stokes equations by the 
mean and fluctuating parts and then time-averaging the equations. This time 
interval must be long compared to the characteristic time of turbulent eddy fluc- 
tuations and short compared to the characteristic time of the, perhaps, unsteady 
main flow. The resulting equations describe the time-dependent flow field but 
are not capable of capturing fluctuating turbulent eddies of any length scale; i.e., 
all turbulent eddies must be modeled. In large eddy simulation, the small scale 
eddies are modeled while the large eddies are computed. Solving the full Navier- 
Stokes equations about a complex geometry at high Reynolds number without 
any averaging is yet beyond the reach of present day and the near future super- 
computers. For this reason, the “Reynolds-averaged Navier-Stokes” equations in 
the “thin-layer” form are used in this project. Therefore, it is necessary to cite 
some of the significant steps taken in this relatively new area. 

Early work started with the computation of laminar flow which is governed 
by the steady Navier-Stokes equations. Although a good survey on laminar flow 
research has been done by Peyret and Viviand [6], certain notable efforts are listed 
below: The shock wave laminar boundary layer interaction by Mac Cormack [7], 
the computation of laminar flow over a compression corner by Carter [8], and 
hypersonic flow over blunt bodies by Jain and Adimurty [9]; also, incompressible 
flow over bluff bodies and airfoils by Thompson et al. [10] and Hodge and Stone 
[11], and supersonic two-dimensional flow over a blunt body with impinging shock 
wave by Tannehill et al. [12]. In three-dimensional laminar flow simulation, we 
can mention the computation of the flow over an inclined body of revolution 
by Li [13], and the laminar flow over three-dimensional compression corners by 
Shang and Hankey [14] and Hung and Mac Cormack [15]. 

Since most aerodynamic flows encountered in actual flight conditions are 
turbulent, the Reynolds-averaged Navier-Stokes formulation draws much of the 
attention from engineers and scientists. Early turbulent flow simulations started 
with two-dimensional flows and include the shock wave/turbulent boundary layer 
interaction study by Wilcox [16] and Baldwin and MacCormack [17], and the high 
Reynolds number transonic flow computation over airfoils by Deiwert [18]. Some 
numerical computations were also carried out for supersonic flow over compres- 
sion corners by Wilcox [19], Shang and Hankey [20], and Hung and MacCormack 
[21] and axisymmetric afterbody computations were carried out by Holst [22]. 
Some of the turbulent simulations yielded very interesting results which closely 
modeled the actual physical phenomena: For instance, the unsteady buffeting 
of a thick circular arc airfoil in transonic flow was simulated by Levy [23], the 
transonic drag polar and lift curve for a supercritical airfoil with buffeting forces 
were reported by Deiwert and Bailey [24], and the simulation of transonic aileron 
buzz was achieved by Steger and Bailey [25]. 
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Three-dimensional solutions of turbulent flows with the Reynolds-averaged 
equations are relatively few. The first simulations were the flows around relatively 
simple geometries; for example the three-dimensional swept shock wave/ turbulent 
boundary layer interaction by Hung and MacCormack [26], and the flow along 
a corner for both transitional and fully turbulent flow by Shang et al. [27]. 
Another notable work was the simulation of a separated flow over a hemisphere- 
cylinder body at an angle of attack by Pulliam and Steger [28]. In most of these 
simulations, the “thin-layer” approximation to the Reynolds-averaged Navier- 
Stokes equations devised by Baldwin and Lomax [29] has been utilized. Very 
recently, computation of more complex geometries involving more sophisticated 
flow physics has been attacked. Work in this area includes the transonic fuselage 
and forebody flows of Cosner [30], the hypersonic wing/fuselage interaction of 
Shang [31], the supersonic blunt-fin/wall interaction of Hung and Kordulla [32], 
the high-subsonic turret flow simulation of Purohit et al. [33], the transonic wing 
flow of Mansour [34], the transonic afterbody flows of Beiwert and Rothmund 
[35] and Deiwert et al. [36], the transonic wing flows of Agarwal and Deese [37], 
the transonic forward-fuselage flow of Chaussee et al. [38], and the high-subsonic 
delta wing and low supersonic, shuttle-like flows of Fujii and Kutler [39]. 

1.2 Motivation 

In the previous section, the role of Computational Fluid Dynamics as a re- 
search and design tool in aeronautical sciences has been introduced. It was also 
stated that there emerged numerous pioneering works in computational Navier- 
Stokes technology. Some of the reasons for this movement are explained as fol- 
lows: First, the computer hardware(memory and speed) have dramatically im- 
proved by the introduction of the so-called supercomputers such as the Cray 1, 
Cray XMP, CDC Cyber 205 [4]. Second, there is an increasing demand from 
the aircraft industry to utilize more realistic approximations in design prob- 
lems. For instance, modern military aircraft must be able to fly in the transonic 
regime at high angles of attack with significant regions of separated flow, and 
the full Navier-Stokes equations are the proper model equations for this pur- 
pose. These reasons have given a tremendous impetus to viscous transonic flow 
research, especially focused on flow separation and shock/boundary-layer inter- 
action. In an excellent survey about transonic Navier-Stokes technology, Mehta 
and Lomax [40] list the flow conditions for which the Navier-Stokes equations ap- 
pear to be required: 1) Shock/boundary layer interaction with no separation, 2) 
shock-induced turbulent separation with immediate reattachment (shock-induced 
separation bubble), and 3) shock-induced separation without reattachment. The 
proper treatment of shock/boundary layer interaction, at least locally, necessi- 
tates the use of the Navier-Stokes equations [41], There are primarily two mo- 
tivations for understanding separated flows [40]. Since uncontrolled separation 
causes stall, controlling and minimizing the effects of separation are desirable. 
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Separation can also be used to improve aerodynamic performance; for instance, 
the rolled up vortex sheet from the sharp leading edge of a delta wing or strake 
creates so-called “nonlinear lift” which can be intelligently exploited by fighter 
aircraft designers. 

Viscous transonic research has been carried out by either viscous-inviscid 
interaction procedures or Navier-Stokes procedures. The viscous-inviscid inter- 
action procedure is quite convenient for the study of separation because the thin 
shear layer approximation remains valid in mildly separated flows [3]. Computa- 
tionally, interaction codes require significantly less computer time than Navier- 
Stokes codes for problems with small regions of separation. However, Navier- 
Stokes formulations are generally much easier to code [42]. Moreover, for strong 
viscous-inviscid interactions, the current boundary-layer based schemes tend to 
break down. The present state-of-the-art associated with three-dimensional in- 
teraction methods applied to separated flows is far from mature [42]. There 
is some progress on computation of turbulent flows over wings primarily using 
integral boundary layer schemes, but these schemes cannot easily compute sep- 
arated flows. The interested reader is directed to Refs. [42], [43], [44] for more 
information. 

The primary difficulties associated with Navier-Stokes procedures are now 
restated for emphasis: 1) computer speed and memory limitations, 2) difficulties 
in turbulence modeling, and 3) difficulties in geometry definition and grid gen- 
eration around sophisticated geometries. The Navier-Stokes simulations cited in 
the previous section utilized relatively coarse grids and required large amounts of 
computer time even on the latest supercomputers. However, the idea of dividing 
the control volume around the geometry into different “blocks” or “zones” and 
solving the governing equations interactively seems very attractive, because this 
technique reduces the difficulties of the first and third items listed above. This 
idea may be elaborated on as follows: Instead of putting the complete data base 
into the main computer memory, only the data associated with one grid zone can 
be put into main memory while the data from the other zones reside on some 
other secondary storage device. This relaxes the computer main memory burden. 
Second, and more importantly, this scheme relieves the task of generating grids 
around geometrically complex domains such as complete aircraft configurations. 
In this fashion, the total flow domain can be divided into multiple zones each of 
which may be in a different form to respond to the geometry, boundary condition, 
and physical constraints. 

In this work, the Euler/Navier-Stokes equations using a zonal grid tech- 
nique are solved to study three-dimensional flow separation on wings. For this 
purpose a computer code called Transonic Navier Stokes (TNS) has been devel- 
oped [45], which presently utilizes four grid zones for isolated wing geometries. 
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In this code, the “thin-layer” approximation to the Reynolds-averaged Navier- 
Stokes equations is solved by utilizing highly clustered grids near the body. Away 
from the body, the inviscid-rotational features of the flow are simulated utilizing 
the Euler equations. To the authors’ knowledge, the TNS computer program 
represents the first three-dimensional Euler/Navier-Stokes zonal algorithm. 
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CHAPTER 2 


DESCRIPTION OF THE ZONAL PHILOSOPHY 


2.1 Introduction 

Historically the zonal philosophy is not new. It was used implicitly in the 
introduction of Boundary Layer theory by Prandtl [3] where the flow was divided 
into two different regions, inviscid and viscous. There are numerous examples of 
viscous-inviscid interaction research which use a type of zonal approach. Also, 
the zonal philosophy finds an application in turbulence modeling. As described 
by Kline et al. [46], zonal turbulence modeling is an attractive alternative to 
universal turbulence models. 

The extensive use and application of zonal schemes in different forms has 
gained popularity recently. Boppe [47] illustrated that using nested grid sys- 
tems around complex geometries such as arbitrary wing-body configurations with 
winglets, pods, canards and tails makes it very practical to construct transonic 
codes. Also, Lee [48], Lee and Ruppert [49], Atta [50] and Atta and Vadyak [51] 
have used the block grid approach for transonic potential computations. Other 
zonal grid approaches for the Euler equations have been used by Benek et al. 
[52], Hessenius and Pulliam [53], Rai [54] and Hessenius and Rai [55]. 

Numerical grid generation in three dimensions is one of the biggest pacing 
items in Computational Fluid Dynamics. Numerically generated curvilinear co- 
ordinate systems are commonly used around arbitrary geometries. Body fitted 
systems are especially helpful for the treatment of boundary conditions. In gen- 
eral, numerical grid generation involves the transformation of a complex physical 
space into a geometrically simple computational domain. When the geometry 
of interest in the physical space is simple enough, the transformation can be 
achieved in a straightforward manner. But as the geometry becomes more com- 
plicated, such as the wing/body/strake/tail combination of a contemporary air- 
craft, the issue of grid generation becomes indeed challenging, and sometimes 
the domain cannot be mapped into a single block. The practical problems en- 
countered in the generation of a single grid about a complex geometry are as 
follows: Often times, trying to generate a single grid results in a skewed grid 
which causes solution inaccuracies. Another difficulty is obtaining the proper 
level of grid refinement in regions of flow gradients, such as a shear layer, a shock 
wave or a vortex. Moreover, trying to configure an acceptable grid via tuning of 
existing grid generation schemes requires substantial amounts of time and effort. 
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One approach to this challenging issue is to divide the domain into a number 
of geometrically simple subdomains, called “blocks” or “zones.” Fig. 1 illustrates 
the idea of zonal grids for a wing/body/tail combination of a fighter airplane. As 
demonstrated here, each component is mapped into a simple rectangular compu- 
tational box. There are a number of ways to structure zonal grids. In general, it 
is quite difficult to classify zonal grid techniques in a clear and concise fashion. In 
Fig. 2, an approach to zonal grid classification is presented. As suggested, zonal 
schemes can be divided into two groups. The first group is composed of schemes 
with overlapping grid zones, and the second group with zones which do not over- 
lap (interface schemes). In the first group (overlap schemes), the overlapping may 
be either “complex,” i.e., the zonal boundaries don’t match each other in any 
way, or “simple” where the zonal boundaries coincide. The Chimera grid scheme 
of Ref. [52] exploits the complex overlap zonal grid technique. The present TNS 
code, however, makes use of the “simple overlap” zonal grid technique. Also, it 
is possible to construct other overlap schemes which blend complex and simple 
overlap schemes. 

The second group of zonal schemes consists of “patched” grids which in- 
terface each other at common boundaries. This group can be divided into two 
subsets according to whether the grid lines at the common surface are continuous 
or discontinuous. In the first scheme, no discontinuity in function, spacing, or 
slope is permitted. The second scheme consists of discontinuous grids. The dis- 
continuity can be in the grid spacing, in the function values or in the slope of the 
grid fines. Also, combinations of different kinds of discontinuities can be used in 
the same grid. This generally occurs when the flow field is divided into different 
zones beforehand, and then the grid in each zone is generated independently (see 
Refs. [48], [49], [56]). The discontinuous grids in Fig. 2 are also “metric discontin- 
uous” in the context of Rai’s [54] definitions. Interface boundary conditions are 
implemented between zones by either interpolation or other special techniques. 
For example, a conservative zonal boundary treatment for discontinuous grids is 
given in Refs. [53], [54], and [55]. 

The advantages of zonal grid schemes are summarized as follows: 

1) The grid for any subregion of the domain can be generated easily. 

2) The grid can be adapted to the problem. This is achieved by placing 
refined zones in regions of high gradients such as shear layers, shock waves, 
jets and wakes, vortices, etc., whereas coarse zones are used in regions of small 
gradients. 

3) Different equation sets can be used in different zones; for instance, the 
Navier-Stokes equations are used in regions where viscous and/or heat transfer 
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effects are important and the Euler equations are used where inviscid rotational 
effects are important. 

4) With the use of grid zoning, an efficient mesh topology can be improved, 
or an inefficient topology can be made acceptable. 

A measure of grid efficiency, the Mesh Efficiency Ratio (MER) was defined in Ref. 
[57], It is the ratio of the number of grid points on the surface of interest to the 
number of grid points in an average two-dimensional surface of the same three- 
dimensional grid. The latter quantity is simply the two-thirds power of the total 
number of grid points. Larger values of MER, in general, indicate more efficient 
grid topologies as opposed to the smaller values of MER which represent less 
efficient topologies. MER values on the order of one for an inviscid grid and on the 
order of one-half for a viscous grid are about average. Of course, this number can 
be changed by up to a factor of two depending on grid clustering. Representative 
values of MER for a few well-known isolated-wing grid topologies are shown in 
Table 1 for both viscous and inviscid grids. This table was constructed with 
conservative stretching estimates for all directions. Thus, larger values of MER 
could be obtained with more rapid stretching. More on the efficiency of mesh 
topologies can be found in Refs. [57], and [58]. 

Table 1. Comparison of several isolated-wing mesh topologies using the MER 
efficiency parameter 


Mesh topology 

MER (inviscid) 

MER (viscous) 

H-H 

0.38 

0.24 

C-H 

0.50 

0.32 

C-0 

0.79 

0.50 

O-H 

0.79 

0.50 

0-0 

1.26 

0.79 

H-H (zonal) 

1.13 

0.71 


As seen from Table 1, it is difficult to advocate the H-H topology about a 
wing involving viscous flow calculations because the boundary layer clustering 
at the wing surface must be continued upstream and downstream of the wing as 
well as outboard of the wing tip. Thus, many points are wasted, and the MER 
is as low as 0.2-0. 3. With a zonal grid approach applied to the inefficient H-H 
topology, much improved MER values can be obtained as seen by the last entry 
in Table 1. However, the complete aircraft configuration is the ultimate aim of 
CFD, not isolated wings! Generally, the Cartesian nature of the H-H topology is 
an appropriate topology for wing/fuselage geometries. This is because, as seen 
in Fig. 1, the H-H grid topology allows each fuselage cross-section to be fitted 
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with an efficient polar-like grid, whereas the wing cross-sectional grid retains the 
Cartesian-like topology. This is a very flexible topology and can be modified 
to handle more complicated configurations. For example, in addition to the 
wing, other lifting surfaces including canards, strakes and horizontal and vertical 
stabilizers can be handled in a straightforward manner. 

2.2 Zonal Grid Generation 

In TNS, grid generation starts with the generation of a global single-zone 
grid which includes the entire flow field. This grid contains no viscous clustering 
and has an H-mesh topology in both the spanwise and chordwise directions. The 
global grid can be generated from either of two approaches: the elliptic solver 
approach of Sorenson and Steger [59], or the parabolic solver approach of Ed- 
wards [60]. They have the capability of generating suitable grids about isolated- 
wing geometries with either free-air or wind-tunnel wall outer boundaries. For 
the most part, the elliptic solver approach of Sorenson, which generates finite- 
difference grids by solving a set of Poisson’s equations, is used in the present 
study. This grid generation method allows control over both grid point spacing 
along, and normal to, the boundaries and angles with which the grid lines inter- 
sect the boundaries. This control is obtained by specification of the boundary 
conditions and the inhomogeneous terms contained within the Poisson s equa- 
tion. Once the boundary definition and grid characterictics have been set, the 
three-dimensional elliptic code solves the governing equations using a successive 
over-relaxation (SOR) procedure. 

Once the base grid is generated, it is divided into zones utilizing a ‘zoning 
algorithm. In the present TNS code an isolated wing grid is configured with four 
zones. The first grid zone (Grid l) is, by and large, the coarse base grid itself 
with a small three-dimensional domain left open into which the rest of the grids 
and the wing are placed. The second grid zone(Grid 2) fills this open domain, 
with a small overlap region in common with Grid 1 and also another small three- 
dimensional domain for the third and fourth grids. Grid 2 is constructed so as to 
contain twice as many grid points in each spatial direction as the original base 
grid. This refinement of Grid 2 relative to the base grid is accomplished via 
cubic-spline interpolation. 

The final two grid zones (Grids 3 and 4) are designed to capture the upper 
and lower surface viscous effects of the wing, respectively. They occupy the space 
left open by the block of points removed from Grid 2, again with a small region 
of overlap included. Grids 3 and 4 are constructed so as to contain the same 
number of points in both the spanwise and chordwise directions as Grid 2. But 
since they are to capture the viscous effects of the wing, the grid points in the 
normal direction are highly clustered. This four-zone grid topology produces an 
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overall grid with a MER of about 0.8. The outer two grid zones are topologically 
represented in the computational domain as cubes with smaller cubes removed 
from the middle, but the third and fourth grids are represented as simple cubes. 

In order to illustrate this zonal grid topology, a series of figures is presented. 
The wing configuration used here is the so-called “WING C” which is a research 
wing and was designed by a cooperative effort between NASA-Ames Research 
Center and Lockheed-Georgia. The planform view of the WING C grid is pre- 
sented in Fig. 3. The wing has an aspect ratio of 2.6, leading edge sweep of 
45 degrees, taper ratio of 0.3, and a twist angle of 8.17 degrees. Fig. 4 shows 
a perspective view of the grid with the symmetry plane highlighted. From this 
view, the wing surface grid resolution including the treatment of the wing tip can 
be seen. A cross-sectional view of the grid in the vicinity of the wing, showing 
details of grid zones 2, 3 and 4 at the symmetry plane, is displayed in Fig. 5. 
Note that the Navier-Stokes grid expands in thickness from the leading edge to- 
ward the trailing edge to better capture the growing boundary layer. In the same 
plane, an expanded view of the WING C leading edge grid detail is presented in 
Fig. 6. The spanwise cross-sectional view of the wing grid showing part of zones 
2, 3 and 4 is given in Fig. 7. The elliptic grid generation scheme was used for 
the grid just presented in Figs. 3-7, and was designed for free-air computations. 

As was mentioned earlier, the elliptic and parabolic grid generators are ca- 
pable of generating grids suitable for use with either free-air or wind-tunnel wall 
boundaries. Fig. 8 shows a typical grid with outer boundary positions specified 
to coincide with the position of the wind tunnel walls from the NASA-Ames High 
Reynolds Number Channel I [61], The grid is plotted in perspective so that de- 
tails on the upper and lower wind tunnel surfaces, the inflow and outflow planes, 
and the wing-symmetry plane are all visible. This grid, which was generated 
by the parabolic grid generation approach [60], becomes the outer coarse grid 
zone (Grid 1). This grid detail near the wing/symmetry plane juncture has been 
removed. It is in this region that grid zones 2, 3 and 4 are located as described 
earlier. The wing geometry used in the wind-tunnel flow simulation of Fig. 8 is 
shown in Fig. 9. It is composed of NACA 0012 cross-sections and has a taper 
ratio of 1.0, a leading-edge sweep of 20 degrees, an aspect ratio of 3, and is rigged 
in the tunnel wall grid at 2 degrees angle of attack. 

The total number of grid points associated with the WING C and NACA 0012 
grids for free air simulations is 165,321. The individual grid point breakdown for 
each zone is as follows: Grid 1, 63 x 26 x 25 = 40,950; Grid 2, 69 x 29 x 21 = 
42,021; Grid 3, 61 x 27 x 25 = 41,175; and Grid 4, 61 x 27 x 25 = 41,175. The 
wind-tunnel wall grid for the NACA 0012 configuration has fewer grid points; it 
is an exact subset of the NACA 0012 wing free-air grid and consists of 149,071 
points (see Fig. 8). 
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2.3 Zonal Interface Scheme 


The zonal approach introduces a new aspect associated with boundary condi- 
tions, namely the interface conditions between neighboring blocks. Whereas the 
interior points of each zone are updated using the standard integration scheme, 
the boundary points are updated via interpolation. It was pointed out earlier 
that the zonal scheme followed here is a “simple overlap” one, and owing to this 
feature, the interpolation procedure is greatly simplified. This is because the 
grid zones are carefully constructed from a base grid such that surfaces requir- 
ing interpolation are coincident. The most complicated zonal interface boundary 
condition involves only a series of one-dimensional linear interpolations. This 
situation is illustrated in Fig. 10. In this hypothetical case, grid Zone 1 is a 
coarse grid which interfaces with Zone 2 which is a fine grid. To facilitate the 
implicitness of the iteration algorithm in the present approach, the blocks are 
overlapped by one to several mesh-cell widths. 

In the interface boundary condit ion scheme, the idea is to obtain flow field in- 
formation from the interior of Zone 1 and to use it to satisfy the proper boundary 
conditions on the ABCD surface of Zone 2. Since the grids of Zone 1 and Zone 
2 have been constructed to coincide with plane ABCD, the process is greatly 
simplified. The ABCD surface is shown in transformed space in Fig. 10b. The 
points from both grid zones are indicated. Note that the grid points from Zone 
1 represent a subset of the Zone 2 points. Hence, the dependent variables at 
the points marked with plus symbols (Zone 2) are computed by a series of linear 
interpolations using the points marked with circular symbols (Zone 1). The pro- 
cess of interpolating information back from Zone 2 to Zone 1, which is required 
for the EFGH surface is even simpler than in the first case. This is because 
the interpolation simply becomes “injection.” That is, even though interpolation 
is carried out, the end result of this operation is the straight transfer of values 
from Zone 2 to Zone 1 without interpolation errors. This is because every grid 
point in the Zone 1 boundary (EFGH), has an exactly matching point on the 
corresponding internal plane (EFGH) of Zone 2. This is readily seen by looking 
at Fig. 10b and noting that every Zone 1 point (circular symbols) has an exact 
matching point in Zone 2 (plus symbols). 

The present interface scheme is organized such that only the two planes 
involved in the interpolation need to be defined: the base and target planes. 
These interface surfaces are referred to as “planes.” However, this is true only in 
the computational domain. In the physical domain these “planes” are actually 
curved three-dimensional surfaces. Each boundary plane at a grid zone interface, 
which requires data for boundary conditions, is called a target plane; and each 
interior plane which supplies interpolation data to a target plane, is called a base 
plane. For instance in Fig. 10, when we want to obtain flow field information 
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from the interior of Zone 1 and to use it as the boundary condition for Zone 
2, the ABCD plane becomes the base plane in Zone 1 and at the same time it 
becomes a target plane for Zone 2. These base and target planes are grouped 
into three categories : J planes, K planes and L planes. Each one is defined by 
six integers. For example, the J planes are defined by : 1) the grid zone number 
from which the base or target surface is obtained, 2) the J plane number in that 
particular zone, 3) the starting K index, 4) the ending K index, 5) the starting L 
index, and 6) the ending L index. For instance in Fig. 10, the ABCD and EFGH 
planes are typical J planes. The K and L planes are defined in a similar way. 

Whereas the base and target plane arrays completely define the dependent 
variable values to be used in the interface interpolation process, the independent 
variable values (spatial coordinates) are obtained from the grid distribution which 
passes through the origin of the coordinate system (the point at the wing root 
leading edge). Obviously, these grid distributions do not uniquely represent all 
the possible grid distributions that may exist in a general three-dimensional mesh. 
However, they offer a very flexible and easy to implement interpolation scheme 
with a reasonably small error. Moreover, if the interpolation is, for example, a 
simple halving of the mesh which is the case for the most part of the TNS grid, 
it is not important to know precise grid coordinates. 

At the beginning of each iteration for grid zone NZ, boundary values for all 
grid zone NZ interface planes are updated from a storage array called BCBUF 
which must be shuffled back and forth from the extended storage. After the 
interface planes and any other standard boundary conditions are updated, the 
interior solution is updated by using the standard integration algorithm. Then 
any base planes that reside in grid zone NZ are interpolated to produce the 
corresponding target planes for use in neighboring grid zone interface boundaries. 
These target planes are stored in the BCBUF array. This completes the iteration 
for grid zone NZ and the algorithm proceeds to grid zone NZ+1. 

A zonal boundary scheme in an ideal sense should be conservative and should 
permit distortion-free movement of flow discontinuities such as shocks and slip 
surfaces across zonal interfaces. Let us investigate our zonal boundary scheme 
with regards to this condition. The upper and lower grid interfaces between the 
inner Euler (zone 2) and Navier-Stokes blocks (zones 3 and 4) are one-to-one 
in both the f and rf directions (streamwise and spanwise directions for a wing 
topology), but in the f direction these planes involve a spacing discontinuity or 
“metric discontinuity.” This feature can be seen most easily in the grid plots 
of the Figs. 4-7. This is important because the interpolation process utilized 
automatically becomes an “injection” for these interface boundaries and therefore 
is conservative. Thus, a strong shock wave crossing these boundaries, which is 
quite likely to occur, should do so with minimal numerical inaccuracies. The 


13 


other interface boundary planes are composed of discontinuous grids and, in 
general, are not conservative, but could be made conservative if necessary by 
using the approach of Rai [54], and Hessenius and Rai [55]. 

2.4 Data Management 

Once the grid is generated and divided into the proper zones, the flow solver 
is initiated. The iteration procedure starts in the outer Euler blocks (first Grid 1, 
then Grid 2) and ends with the two Navier-Stokes blocks, first the upper Navier- 
Stokes block (Grid 3) and then the lower block (Grid 4). Only one iteration 
using a spatially varying time step is completed in each zone before passing to 
the next. However, many iteration and/or time-stepping strategies could be used 
to improve convergence. Only the flow field solution (Q arrays), transformation 
Jacobian (J), metric quantities, and the turbulence model arrays (when appro- 
priate) associated with a single block reside in main memory of the Cray-XMP 
at a time. The information associated with the other blocks resides in extended 
storage. On the Cray XMP this device is called the Solid State Device or SSD. 
The SSD is utilized functionally in the same manner as standard rotating disk 
extended storage. However, the SSD extended storage is physically made up of 
semiconductor memory and therefore is much faster. Using the SSD instead of 
disk greatly reduces I/O wait time, and for jobs which are normally I/O bound 
this is a significant advantage. 

In the present zonal approach, the use of the SSD allows a great deal of 
flexibility because a larger number of grid blocks can be supported without sig- 
nificant additions to main memory. Of course the limiting factor with regard to 
this point is the size of SSD. The SSD installed with the NASA Ames Cray XMP 
presently has 16 million 64-bit words of memory. This can easily be extended 
to 32 million words if half precision (32 bits) is used. The current version of the 
TNS code, with grid dimensions as outlined above in the grid generation section, 
requires 5.8 million words of SSD storage. All arrays on the SSD are stored in 
64-bit precision with the exception of the metric arrays which are stored with 
32-bit precision. A test was performed with 64-bit metric storage that produced 
results very close to the case with 32-bit metric storage. 

To allow more space in main memory, the metrics are shuffled into main 
memory from the SSD in two-dimensional planes as needed. This allows the 
maximum grid size of each grid zone to be about 50,000 points. Because the flow 
solver algorithm used in TNS is an Alternating Direction Implicit (ADI) algo- 
rithm with implicit sweeps in all three directions, the metrics must be transferred 
into main memory with three different orientations, in x-y planes, x-z planes and 
y-z planes. Thus, there are three different metric orientations stored on the SSD. 
Because of the availability of so much storage on the SSD, this does not cause 
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any problem and makes overall memory management more efficient. In addition, 
each of the metric arrays is required in main memory several times for each grid 
zone during each iteration. This places extreme demands on the TNS I/O re- 
quirements. Nevertheless, with the efficiency level afforded by the SSD, these 
I/O requirements are handled without problem. 
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Fig. 1. Body-conforming zonal grids of a fighter aircraft configuration in the physical 
and computational spaces. 
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Fig. 2. An approach to the classification of the present zonal grid techniques. 
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Fig. 4. Perspective view of embedded grid with upper symmetry plane (y = 0, z > 
and wing surface highlighted. 


19 



GRID 



20 



X 


Fig. 6. Expanded view of the WING C leading edge grid details. 
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Fig. 7. Spanwise cross-sectional view of the wing grid showing parts of Zones 2, 3 
and 4. 
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Fig. 8. Global finite-difference grid showing detail on wind tunnel walls. 
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Fig. 9. Plan form view of NACA 0012 wing: AR - 3, \ LE = 20°, a T wi^r = 0°, 
TR - 1 . 
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(a) TWO-ZONE GRID SHOWING OVERLAP AT A BCD AND 
EFGH PLANES IN PHYSICAL SPACE 
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(b) GRID POINT DETAIL IN THE OVERLAP REGION IN 
TRANSFORMED SPACE. 


Fig. 10. Interface procedure of the zonal grids. 
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CHAPTER 3 


DESCRIPTION OF THE PHYSICAL EQUATIONS 
AND THE NUMERICAL METHODS 


3.1 Introduction 

The equations solved in this study are the Euler and Reynolds-averaged 
Navier-Stokes equations written in strong conservation form. The Reynolds- 
averaged equations are simplified by using the standard thin-layer approximation 
for the viscous terms. Once the zonal scheme is developed, any of the standard 
Navier-Stokes integration algorithms can be used. In this study, the basic gov- 
erning equations and numerical algorithm including the turbulence model have 
been taken from the Pulliam-Steger ARC3D computer code [62], To establish a 
basic level of understanding for the present approach, certain characteristics of 
the ARC3D computer code used within TNS are now briefly discussed. 

The governing equations are generally nondimensionalized by free-stream 
quantities and are transformed to the computational domain (£,ri,$) so as to 
preserve the strong conservation form of the equations (see Sec. 3.2). The time- 
dependent metrics which occur in the most general form of the Euler/Navier- 
Stokes equations have been omitted. This was done to save computer memory. 
Thus, the present formulation is incapable of time varying grids but still has the 
capability of resolving time-dependent flow-field physics. However, adding the 
time-dependent grid capability back into the present code is not a difficult task. 

Two numerical algorithms have been investigated within the TNS computer 
code, an ADI algorithm which solves block-tridiagonal matrices along each coor- 
dinate direction and a diagonalized algorithm which solves scalar pentadiagonal 
matrices along each coordinate direction (see Sec. 3.3). The first algorithm is 
a variation of the Beam- Warming ADI scheme (63], and the second is an ex- 
tension of this scheme due to Pulliam and Chaussee (64). Both schemes use the 
standard second-order-accurate central differencing of the governing equations to 
construct the appropriate spatial differencing scheme. The block scheme uses a 
fourth-order-accurate smoothing operator on the right-hand side of the iteration 
algorithm and a second-order-accurate smoothing operator on the left-hand side. 
The diagonal scheme uses a fourth-order-accurate smoothing operator on both 
the left- and right-hand sides. 

The diagonal algorithm has been implemented with two spatially varying 
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time-step philosophies; one which scales the time step with the transformation 
Jacobian and the other with a combination Jacobian/local-solution variation. 
The block algorithm uses a time-step which is constant in a spatial sense as well 
as from grid zone to grid zone. 

All of the results presented in this report have been computed using the di- 
agonal algorithm with fourth-order implicit smoothing and the Jacobian-scaled 
time-step variation. This algorithm combination seemed to be the most compu- 
tationally efficient of all the variations tested. More discussion of these aspects 
can be found in Flores [65] where a comparison of these algorithm variations 
in conjunction with zonal grid variations has been made. The turbulence model 
used in the TNS computer program was the Baldwin-Lomax algebraic model [29] 
because of the ease with which it can be implemented. 


3.2 Governing Equations and Approximations 

3.2.1 Equations in Non-Dimensional Form 

The strong conservation law form of the Navier-Stokes equations are used 
for shock capturing purposes. The equations in Cartesian coordinates in non- 
dimensional form can be written as (cf. Peyret and Viviand [66]) 
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The Cartesian velocity components u, v, and w are nondimensionalized by 
floo (the free-stream speed of sound), density p is nondimensionalized by p ^ ; 
and the total energy per unit volume e is nondimensionalized by p^a^. Pressure 
can be found from the ideal gas law as 

P= {l - 1) (e - 0.5p(u 2 + v 2 + w 2 )} (3.4) 

and throughout 7 is the ratio of the specific heats. Also, k is the coefficient 
of thermal conductivity, p is the dynamic viscosity, and A from the Stokes’ 
hypothesis is —2/3 p. The Reynolds number is Re and the Prandtl number is 
Pr. 


To enhance numerical accuracy and efficiency and to handle boundary con- 
ditions more easily, the governing equations are transformed from the Cartesian 
coordinates to general curvilinear coordinates (Fig. 11) where 
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The resulting transformed equations are not much more complicated than the 
original Cartesian set and can be written in nondimensional form as 
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where [/, V, and H 7 are contravariant velocity components written without metric 
normalization. The viscous flux terms are given by 
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where the components of the shear-stress tensor and heat-flux vector in nondi- 
mensional form were given in (3.3). Here, the Cartesian derivatives are expanded 
in y, f space via chain-rule relations such as 
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Finally, the metric terms are obtained from chain-rule expansion of x^, y^, etc., 
and solved for f x , f y , etc., to give 
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3.2.2 Thin-Layer Approximation 


In high Reynolds number flows, the viscous effects are confined to a thin layer 
near rigid boundaries. In most cases, there are only enough grid points to resolve 
the gradients normal to the body by clustering the grid in the normal direction, 
and resolution along the body is similar to what is needed in inviscid flow. As a 
result, even though the full derivatives are retained in the equations, the gradients 
along the body are not resolved unless the streamwise and circumferential grid 
spacings are sufficiently small. Hence, for many Navier-Stokes computations, 
the viscous derivatives along the body are dropped. This leads to the thin-layer 
Navier-Stokes equations. 

The thin-layer model requires a boundary layer type coordinate system. In 
our case, the £ and r\ directions are along the body and the viscous derivatives 
associated with these directions are dropped, whereas the terms in f are retained 
and the body surface is mapped onto a constant f surface. Thus, Eq. (3.6) 
simplifies to 
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It should be emphasized that the thin-layer approximation is valid only for high 
Reynolds number flows. Also, very large turbulent eddy viscosities invalidate the 
model. 

3.3 Numerical Method 

The finite-difference schemes used are the implicit approximate factorization 
algorithm in delta form by Beam and Warming [63] and a diagonal implicit 
algorithm by Pulliam and Chaussee [64]. Implicit methods are chosen to avoid 
the time-step restriction when small grid spacing is used. When spatial resolution 
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of large gradients such as that associated with shock waves or viscous effects are 
needed, highly refined grids are required. These requirements yield stiff problems 
which in turn restrict the time step to very small values. Implicit methods are 
useful in avoiding such stiffness, and large time steps compared to those of explicit 
schemes can be used Without degrading accuracy. 

3.3.1 Beam- Warming Block ADI Algorithm 

The basic Beam-Warming algorithm is first- or second-order accurate in 
time and second- or fourth-order accurate in space. The equations are factored 
(spatially split) which, for a given time iteration, reduces the process to three 
one-dimensional problems. Due to the second-order central-difference operators 
employed, the algorithm produces block tridiagonal systems for each spatial di- 
mension. The stability and accuracy of the numerical algorithm is described by 
Beam and Warming [63]. According to the linear analysis, the numerical scheme 
is unconditionally stable in two dimensions but in actual practice time step lim- 
its are encountered because of the nonlinear nature of the equations. However, 
this limitation is much less stringent than comparable explicit schemes. In three 
dimensions the algorithm is unconditionally unstable, but stability is maintained 
by the addition of artificial dissipation terms. 

The finite-difference algorithm due to Beam and Warming applied to (3.11) 
results in the following approximate factorization: 

{I + hS^A n + D[ 2) ){I + h6r,B n + Z?J 2) )x 

(/ + M ? C n - hRe~ 1 6sJ~ i M n J + £>j 2) )AQ n = R n ( 3 - 13 ) 

= -A r{6cE n + 6 v F n + 6^G n - Re- l 6^S n ) - D^Q n 

where 6 is the central-difference operator and A and V are forward and backward- 
difference operators, e.g. 

6*Q = | Q{t + Ae,i7,f) - QU ~ Af , 77 , ?)]/2Af 

= |Q(£ + A e,7?A) - gtf,i7,?)]/A£ (3.14) 

= [Q(£,»?,f) - QU - A£,t?,<;)]/A£ 

Indices denoting spatial location are suppressed and h = A r corresponds to first- 
order time-accurate Euler Implicit and h = At/2 to second-order time-accurate 
Trapezoidal Rule. D ^\ and D ^ are the implicit and is the explicit 

smoothing operator which will be explained later, in Section (3.4). 

The Jacobian matrices A n , B n and C n are obtained by linearizing the flux 
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vectors E n , F n and G n in time such that 


E n+ 1 = E n + A n (Q n+1 - Q n ) + 0(Ar 2 ) 

F n+ 1 = F n + B n (Q n+1 - Q n ) + 0(Ar 2 ) 
G n+1 = G n + C n {Q n+1 - Q n ) + 0(Ar 2 ) 
where indices denoting spatial location are suppressed again and 
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are the flux jacobian matrices. These flux jacobians and the viscous coefficient 
matrix M, which comes from the time linearization of the viscous vector 5 n+1 , 
are documented in Appendix A. 


3.3.2 Pulliam-Chaussee Diagonal ADI Algorithm 

Block tridiagonal-matrix inversions constitute the major portion of numeri- 
cal work associated with the standard Beam- Warming algorithm. Equations (3.6) 
are a coupled set of 5 equations and thereby produce a (5 x 5) block-tridiagonal 
structure for the implicit operators of Eqs. (3.13). The diagonal version of the 
standard algorithm due to Pulliam and Chaussee [64] overcomes this difficulty. 
In this algorithm, rather than inverting block-tridiagonal matrices in each direc- 
tion, scalar pentadiagonal matrices are inverted. This is computationally more 
efficient. 


The Jacobian matrices, A, B and C have a set of eigenvalues and a complete 
set of distinct eigenvectors. Similarity transformations can be used to diagonalize 
A, B and C, 

A = T^T^, B = T ri k r ,T- 1 , C = T,k,T~ 1 . (3.17) 
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The similarity transformation matrices T$, T f and their inverse matrices 
are given in Appendix B. Relations exist between T £, and T ^ of the form 

N — Tr'Tr,, = T~ l Tc, P = T ~% , ,P _1 = T~ l T v (3.20) 
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After applying the similarity transformations of (3.17) and identities (3.20) in 
Eq.(3.13) and exchanging the smoothing operators with new ones, the diagonal 
form of the standard algorithm reads 

Ti(I + /i^A e - hDi\z)N{ I + h6r,k n - hD { |„)x 
P{I + h6sks-hDi\s)Tf l £Q n = R n 

The spatial accuracy of the standard and diagonalized algorithms for steady-state 
problems (i.e. AQ n -> 0 as n -+ oo) is determined by the type of differencing in 
forming R n . Since the modification that produces the diagonal algorithm does 
not effect R n , both schemes will have the same steady-state solution assuming 
that the steady-state solution is independent of the convergence path; i.e., that 
the steady state is unique. For constant coefficient matrices A, B and C, the 
diagonal algorithm reduces to the standard algorithm because the eigenvector 
matrices are also constant. Therefore, the linear stability analysis of Beam and 
Warming (63] also holds for the diagonal algorithm. 


By numerical experimentation with the TNS code using both schemes, it was 
found that the diagonal algorithm performs much better in terms of convergence 
rates (see Flores [65]). To further enhance the convergence rate, two spatially 
varying A r schemes have been used. 


A r 


A t re j 

^max 


and 


At 


At re/ 

i + Vj 


(3.23) 
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where J is the Jacobian of the coordinate transformation, \ max is the maximum 
eigenvalue of the flux jacobians A, B and C, and A t re f is a reference time step. 
The latter time stepping scheme proved to be better and was used throughout 
this study. 

3.4 Boundary Conditions 

Since the finite-difference scheme used in the TNS code is implicit, it might 
be expected that the boundary conditions would also be implicit. But the imple- 
mentation of the implicit boundary conditions necessitates a modification of the 
inversion matrices which is not very practical. Moreover, past experience is such 
that explicit boundary conditions have been successfully employed in conjunc- 
tion with the implicit codes. Also, explicit treatment of the boundaries leads to 
a simple and flexible scheme. Boundary conditions become a modular element 
that can be put in or pulled out of a computer program without touching the 
heart of the implicit code. 

Unknown values of Q on the boundaries are updated explicitly, and this 
leads to a zeroth-order time extrapolation. The space extrapolations are either 
zeroth order or first order according to the physical and geometric constraints. 
Since in the TNS code a zonal approach is employed, there are two types of 
boundaries at which conditions must be specified: 1) the physical boundaries 
such as inflow, outflow, and body surface, and, 2) the zonal boundaries where 
“zones” are patched together. Treatment of zonal boundaries was presented in 
Chapter 2. 

At the farfield boundaries, free-stream values are specified. This is shown 
in Fig. 12. At the upper boundary DC, lower boundary AB, side boundary EF, 
and inflow boundary AD, the conditions are free-stream, i.e., Q — Q TO . The con- 
ditions at the outflow boundary are found by a zeroth-order extrapolation from 
the last plane such that Qjmax — Qjmax- i- The symmetry plane bound- 
ary conditions are more complicated: Here a zeroth-order space extrapolation is 
used for the density. A first-order extrapolation is employed for the x-direction 
Cartesian velocity component u and z-direction Cartesian velocity component 
w , while the spanwise y-direction Cartesian velocity component v is set equal 
to zero to force symmetry. Pressure is also found by a first-order extrapolation 
and the total energy is computed using Eq. (3.4). At solid surfaces, the no-slip 
boundary condition is employed for the three velocity components. Surface pres- 
sure is found by using the boundary layer approximation dpjdn « 0 such that 
Pi = P 2 , where p\ is the static pressure on the body surface, and p 2 is the static 
pressure on the first computational surface off the body. The surface density is 
found from an adiabatic wall assumption which yields pi = p 2 . Total energy is 
computed using Eq. (3.4). 
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3.5 Artificial Dissipation Models 

The numerical solution of the Euler/Navier-Stokes equations for flows involv- 
ing captured shock waves necessitates the use of numerical or artificial dissipation 
either by using a scheme which is inherently dissipative or by using a nondissipa- 
tive scheme with specially designed dissipative terms. The former class of schemes 
employ some form of upwind differencing under the assumptions of characteristic 
theory and wave propagation. The schemes due to Steger and Warming [67], Roe 
|68], Van Leer [69], Osher and Chakravarty [70] and Harten’s TVD methods [71] 
all fall in this category. In this class of schemes, the dissipation is automatically 
included without user control. The basic ARC3D algorithm [62] used in the TNS 
code is of the second variety and consists of a centrally-differenced nondissipative 
scheme with user-controlled smoothing terms. The artificial dissipation is added 
for two reasons; first to control the odd-even decoupling of grid points typical of 
central differencing, and second to control strong nonlinear effects such as shocks. 

It can be shown that upwind schemes are equivalent to central difference 
schemes with added dissipation. Pulliam[72] demonstrates that adding a fourth- 
difference dissipation to a central-difference scheme produces the equivalent of a 
second-order upwind scheme, whereas the use of a second-difference dissipation 
term produces a first-order upwind equivalent. There are different dissipation 
models utilized in the block and diagonal algorithms described earlier in Sec. 
(3.3). 

The block algorithm which was given in Eq. (3.12) uses a constant coefficient 
dissipation model. In this approach, an explicit fourth-difference dissipation 

D {4) Q n = e e Ar J -1 [(V ^A^) 2 + (V„A„) 2 + (V, A ? ) 2 ]JQ n (3.24) 

is added to the right-hand side of the equation and an implicit second-difference 
dissipation 

d[ 2) = -e.-J-WVfcJ, * = (3.25) 

is inserted into the respective implicit block operators. Based on linear stabil- 
ity theory, the use of the explicit dissipation alone produces an explicit stability 
bound. If the second-difference implicit dissipation is added, the algorithm be- 
comes unconditionally stable if is sufficiently large (cf. Ref. [28]). The proce- 
dure is to set c e = A t and e, = 2e e which results in a consistent definition of e e 
such that as the time step increases, the amount of artificial dissipation added 
relative to the spatial derivatives of convection and diffusion remains constant. 
This is because all of the convection and diffusion terms in the approximate fac- 
torization algorithm of the Eq. (3.13) are factored by the time step A t. Instead 
of second-difference implicit dissipation, using fourth-difference implicit dissi- 
pation would improve convergence. However, this would require the inversion 
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of extremely expensive block pentadiagonal matrices rather than the relatively 
simple block tridiagonal matrices [72], [73]. 

The diagonal version of the implicit approximate factorization algorithm 
due to Pulliam and Chaussee, which was described in Sec. (3.3), substantially 
reduces this computational cost. The more appropriate fourth-difference dissipa- 
tion terms can be added, and only require the inversion of scalar pentadiagonal 
matrices as opposed to block pentadiagonal matrices in the block algorithm. 
The diagonal algorithm was given by Eq. (3.22). The combination nonlinear 
dissipation operators are included as follows: 


-DtU - ^ _ Z&tlJ (3.26) 


( 2 ) 


.(4) 


The terms for Di\ n and D t | ? have a similar form. The same dissipation model 
has been implemented in the right-hand explicit side of the algorithm. Hence the 
constant-coefficient fourth-difference dissipation operator in Eq. (3.13) is 

replaced by a combination nonlinear dissipation operator 


D e \‘ (3.27) 


The dissipation operators defined by Eq. (3.24) and (3.25) include the metric 
Jacobian J so that they act upon the unsealed Q and not on Q. Hence the 
dissipation is not directly influenced by the grid distribution because smooth 
functions of Q may occur in regions of rapid mesh variation. This form of 
smoothing operator is inspired by the recent work applying flux limiters to up- 
wind schemes and TVD concepts which suggest that the best approach for an 
upwind algorithm is to use locally first-order upwind differences at the shock 
and second-order elsewhere[72]. This is accomplished by using switch operators 
cjffc anc * e } 4 ,k which govern the type of smoothing, e.g. second-difference smooth- 
ing near shocks and fourth-difference smoothing elsewhere. The term o }t k is a 
spectral radius scaling which in three dimensions is defined as 


o it k,i = \U\ - a^kl + Jl + Jl 

+ \v\ + a \fnl + + nl (3.28) 

+ |W| + « v te + tf + ? 2 

and is a sum of the spectral radii of the flux jacobians A, B and C. The interested 
reader is referred to Pulliam [72] for a good review of this subject. The effect of 
different smoothing models in the actual computations for the TNS code will be 
presented in Section (4.3). 
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3.6 Turbulence Model 


Solving turbulent flows numerically is a very difficult task due to the extreme 
time and space scales associated with turbulent motion. Hence, the main method 
of approach in Computational Fluid Mechanics and Heat Transfer for solving 
turbulent flows is through the Reynolds averaged Navier-Stokes equations. Time 
or “Reynolds” averaging the equations of motion give rise to new terms which can 
be interpreted as “apparent” stress gradients and heat flux quantities associated 
with the turbulent motion. The apparent turbulent stresses in compressible flow 
can be written as [74] 

{Tij)turb = -puty ( 3 - 29 ) 

and the apparent turbulent heat flux components as 

Q 

“(V .q)t«rfc = -^{pc p T^) (3.30) 


In order to predict turbulent flows by applying finite-difference methods to 
the Reynolds equations, it is necessary to make assumptions for the terms in 
Eqs. (3.29) and (3.30). Boussinesq (1877) suggested that the apparent shearing 
stresses might be related to the rate of mean strain through an apparent turbulent 
or “eddy” viscosity given in incompressible flow by 


- pu[u 
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PT 


dui 

dx-j 


+ 


dun 
dx{ / 


(3.31) 


Closure for the Reynolds heat flux-term, pc r T'uj is treated in algebraic models 
by a form of Reynolds analogy. The Reynolds analogy is based on the similar- 
ity between the transport of heat and momentum and applied to the apparent 
turbulent conductivity in the assumed Boussinesq form 


pCpT'u'j = ~k T ~Q^r (3-32) 

Experiments reveal that the ratio of the diffusivities for the turbulent transport 
of heat and momentum which is called the turbulent Prandtl number, Pr t = 
PT c p/^T ■> I s a well-behaved function across the flow and in algebraic models is 
generally taken to be 0.9. Using the turbulent Prandtl number, the turbulent 
heat flux is related to the turbulent viscosity and mean flow variables as 


-pCpT'u'j 


CpPr &T 
Ptt dxj 


(3.33) 
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Hence, in large scale Navier-Stokes computations an accurate and compu- 
tationally efficient method of computing eddy coefficients is needed. Since the 
algebraic eddy viscosity models proved to be fairly accurate and computationally 
efficient, the Baldwin-Lomax algebraic eddy viscosity model [29] is utilized here. 
In the Baldwin-Lomax model the effects of turbulence are simulated by replac- 
ing the molecular coefficient of viscosity mm with the effective viscosity mm + MT 
in the stress terms of the laminar Navier-Stokes equations. In heat flux terms 
k/c p = p/Pr is replaced by p/Pr t Mr/^fr- 


The Baldwin-Lomax model is a two-layer algebraic model in which Mr is 
given by 


Mr 


\ (Mr) inner? V ^ V crossover 
l (Mr) outer ? V ^ y crossover 


(3.34) 


where y is the normal distance from the wall and y cr ossover is the smallest value of 
y at which values from the inner and outer formulas are equal. The eddy viscosity 
coefficient in the inner layer is based on the Prandtl mixing-length theory 


{Pt) inner — P ^ |w 


(3.35) 


The parameter t is the mixing length corrected with the Van Driest damping 
factor to account for the laminar sublayer 


l - ky [1 - exp(-y + /A+)} (3.36) 

where k -- 0.4, A + = 26, and |w | is the magnitude of the vorticity given by 

M = \/ (“y - Ox) 2 + (o 2 - Wy) 2 + (yj x - u z ) 2 (3.37) 


and 
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-i- p w u T y 

Pw 


\J Pwfwy 

Pw 


(3.38) 


The eddy viscosity coefficient in the outer layer is based on the distribution of 
vorticity which is used to determine the length scale and is given by 


(mt) Outer ~~ K Ccp p FwakeFkleb{ y) (3.39) 


where K — 0.0168 is the Clauser constant and Ccp = 1.6 is an additional 
constant. Fwake is found via 


Fwake = min 


ymaxFmax 


or 


F'WK ymax u \)ipl F n 


(3.40) 
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where Cwk = 0.25 and 


UDIF - (v « 2 + v 2 + w 2 )mai - (\/u 2 + V 2 + tU 2 ) mtn ( 3 - 41 ) 

In Eq. (3.41) the second term is taken to be zero (except in wakes). The quan- 
tities y m ai and F ma x are determined from the function 

F(y) = y[w|[l - exp(~y + /A + )} (3.42) 


In wakes the exponential term is set equal to zero. The quantity Fmax is the 
maximum value of F(y) that occurs in the profile and Umax is the value of y at 
which Fmax occurs. The function Fkleb is the Klebanoff intermittency factor 
given by 


Fkleb{v) — 


1 + 

' yM AX ' 


-1 


(3.43) 


and C kleb — 0.3. 


Although not seen explicitly, the Reynolds number enters into the compu- 
tation of eddy viscosity through the computation of y + . When the variables in 
Eq. (3.38) are nondimensionalized, the following expression is obtained: 


y 


+ 
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Rz PxjjTx i 

Muj 


-y 


(3.44) 


In high Reynolds number flows, as Re -* oo, the length scales / and F(y) in the 
inner and outer layers go to zero, which result in vanishing values of /x j 1 . For 
more discussion on the Baldwin-Lomax turbulence model, the reader is referred 
to the original paper [29]. 


3.7 Computational Aspects of the Zonal Algorithm 

So far the building blocks of the TNS program were presented: namely, 
the zonal algorithm, related data management and the numerical methods and 
flow models. It is very important to explore the characteristics of this code if 
it is intended to be used in routine flow simulations. The speed, convergence 
characteristics, efficiency and reliability of the code must be assessed. 

The convergence rate of the diagonal version of the code (see Sec. 3.3.2) 
versus the block (standard) version (see Sec. 3.3.1) is shown in Fig. 13 (see Ref. 
[65]). The geometry used in this simulation was a NACA 0012 wing (see Fig. 9). 
The free-stream flow conditions were M 0 0 — 0.826, a = 2° with a Reynolds 
number based on chord of 8 million. This is a moderately difficult case involving 
a strong shock with a moderate amount of separation. In Fig. 13, the abscissa 
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denotes the number of iterations for the algorithm, and the ordinate denotes the 
L2-norm of the residual. The time step used in the block version was At = 0.004, 
where a unit time interval denotes the time required by a particle moving at the 
free-stream speed of sound to travel one chord length. This was the largest 
time-step possible while still maintaining stability of the code. 

The diagonal version used a variable time-step (see Eq. 3.23) and thus 
received convergence accelereration from this aspect. However, despite this, the 
main speed-up associated with the diagonal algorithm is not due to the variable 
time-stepping alone, but the variable time-step coupled with the correct implicit 
treatment of the numerical dissipation terms. Pulliam and Steger [73] reported 
that the variable time-step had often worked poorly in viscous flow until the 
second-order implicit dissipation terms were replaced by the fourth-order terms 
which are consistent with the fourth-order explicit dissipation of the right-hand 
side of the algorithm (see Eqs. 3.26 and 3.27). The same replacement is, of 
course, possible in the block algorithm, but would necessitate the inversion of 
expensive block pentadiagonal matrices. 

The slow rate of convergence associated with the block scheme seems to occur 
in the outer zones. The viscous-zone block-scheme residual drops rapidly during 
the first 1000 iterations but then flattens out. At 5000 iterations, the residuals 
associated with all zones have dropped by about one to two orders of magnitude. 
In contrast, the convergence rate of the diagonal scheme drops rapidly in all 
zones. A three-order of magnitude drop in the L2-norm for a 150,000 grid point 
mesh occurs in 500 iterations which costs 2450 seconds of CPU time. This is 
equivalent to 4.9 seconds of CPU time per iteration and 3 x 10 -5 seconds of 
CPU time per grid point per iteration. The convergence rate of the diagonal 
algorithm is nearly 40 times that of the original block algorithm. This is due to 
the proper linearization of the dissipation terms and variable time stepping as 
stated earlier, coupled with the decrease in arithmetic operation count owing to 
the diagonal algorithm. 

The convergence characteristics of the algorithms as discussed above and 
other computational aspects of the TNS code were thoroughly explored by Flores 
[65], It was concluded that: 1) although convergence rates are vastly different 
for different algorithms, the solution at steady state is not affected, 2) the use 
of explicit boundary conditions at zonal boundaries has no detrimental effect 
on the convergence characteristics of the algorithms, 3) there was no noticeable 
improvement in the convergence rate as a result of changing the extent of overlap 
between zones, and 4) use of different CFL values in different zones increases the 
convergence rate by 10-20%. 

Data management of large application codes is an important issue which 
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must be handled carefully. Memory, speed and I/0(input/output) requirements 
of the TNS code are quite substantial. The key element for handling the memory 
and I/O requirements of the TNS code is the SSD (Solid State Device) of the 
NASA Ames Cray X-MP (see Sec. 2.4). The SSD has 16 million 64-bit words 
of memory. It will be instructive to give some statistics about how this code 
performs in actual situations, and in particular, how the data management of 
the metrics was achieved. 

To exploit the ADI character of the flow solver algorithm, which necessitates 
implicit sweeps in all three directions, the metrics are transferred from SSD into 
main memory with three different orientations: x—y planes, x — z planes, and y — z 
planes. Hence, there are three different copies of the metrics stored on the SSD 
corresponding to the three different orientations, which is not a problem because 
of the ample storage space in the SSD. However, the number of grid points for 
an adequate grid resolution is still limited by the capacity of the main memory 
(1 million 64-bit words on Cray X-MP). Each of the metric arrays is required in 
main memory several times for each grid zone during each iteration. This places 
extreme demands on I/O requirements. Nevertheless, the SSD handles these 
I/O requirements without a problem. The computational statistics displayed in 
Table 2 demonstrate this for two cases (Moo — 0.80 and 0.95 and Re = 8 x 10 7 , 
a = 5°). The grid used for both calculations is similar to that presented in 
Fig. 8 and consists of four zones with a total of 166,621 grid points. As seen 
from Table 2, the higher Mach number case converges more slowly than the 
lower Mach number case. This is not surprising, since the Mqo = 0.80 case is 
only moderately separated and the Moo = 0.95 case is massively separated. To 
advance a typical interior grid point in the TNS program one time level requires 
approximately 2030 floating point operations. Thus, for the cases displayed in 
Table 2, the TNS program execution rate on the Cray X-MP computer is about 
63 MFLOPS (million floating point operations per second). 

From the statistics in Table 2, the efficiency level associated with the SSD 
I/O is quantitatively established. During the Moo = 0.95 case, the program 
transfers almost 24 billion 64-bit words between the SSD and main memory 
using about a quarter of a million read/write requests. The I/O time charge for 
this case is just over two minutes. The estimated I/O time charge for this case, 
assuming disk is used to replace the SSD, is over 24 hours. (Because the I/O 
utilized in the TNS program is asynchronous, these I/O time charges represent 
only that portion of the I/O time that was not “covered up” by the CPU during 
execution.) For a detailed survey of the computational aspects of the TNS code, 
please refer to Ref. [75], 
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Table 2. Computational Statistics from the TNS Program 

(NACA 0012 wing, Re c = 8 x 10 7 , AR = 3.0, TR = 1.0, A = 20°) 


Quantity 

Moo = 0.80 

Moo = 0.95 

Iterations* 

825 

1890 

CPU Time (Hrs) 

1.25 

2.82 

Read/ Write Requests (Memory /SSD) 

110550 

253260 

64-Bit Words Transferred (MWDS) 

10391 

23806 

SSD I/O Time (Sec) 

61.5 

140.8 


Estimated Disk I/O Time (Hrs) 10.7 24.4 

(If disk had been used in place of SSD) 


* Avg residual reduced three orders of magnitude (all grid zones). 
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Fig. 12. Specification of the boundary conditions. 
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CHAPTER 4 


SIMULATION OF WING FLOW FIELDS 


4.1 Introduction 

Understanding the nature of separated flows is important for an aircraft 
designer. If the aerodynamic design of a flight vehicle with flow separation is to be 
successful over the full range of flight conditions, the vehicle must be controllable 
at all times and possess no rapid changes in force and moment characteristics. 
Lifting surfaces such as wings, canards, strakes, etc., play an important role in 
the design process. Accounting for separated flows about these surfaces is one 
of the most crucial tasks of the aerodynamic design. Designers either want to 
reduce the adverse effects of separation, or exploit separation for useful purposes. 
The latter, for instance, is the case for a rolled-up vortex sheet from the sharp 
leading edge of a delta wing or a strake of a fighter aircraft, in which case the 
vortex provides so-called “nonlinear lift.” 

Our present understanding of three-dimensional separation about wings 
comes principally from flow visualization experiments with surface-oil techniques 
in wind tunnels, smoke tunnels, and dye flow in water tunnels. However, there 
have also been some attempts to compute separated flows on wings in recent 
years. Numerical simulations of the leading edge vortex phenomenon associated 
with a delta wing have been made using potential flow procedures such as panel 
methods or vortex lattice methods [76], [77]. Other researchers have used the 
Euler equations to solve the same problem [78], [79], [80], [81]. For a discussion 
about the prediction of separation, which is a viscous phenomenon, using the 
Euler equations, see Ref. [82]. 

Another technique used in the computation of separated flows is based on 
viscous-inviscid interaction methods (see Sec. 1.2). These methods are quite 
appropriate for computing mild separations but, at present, tend to break down 
for large streamwise separations. Some progress has been made on computing 
turbulent flows over wings primarily using integral boundary layer schemes, but 
these schemes cannot easily compute separated flows [42], [43], [44]. 

Because of the shortcomings of the methods described in the previous para- 
graph, the Navier-Stokes equations seem to be the proper model equations for 
computing massively separated flows on wings [40]. This line of research is quite 
recent, with only a few publications on the subject. The computation of the 
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supersonic flow over a blunt delta wing by Vigneron et al. [83], the leading-edge 
separation vortex over a delta wing at high angle of attack by Fujii and Kutler 
[39], and the simulation of a tip vortex off a low-aspect-ratio wing at a transonic 
speed by Mansour [34] are a few examples. Also very recently, transonic wing 
solutions have emerged due to Agarwal and Deese f 37], Vadyak [84] and Holst et 
al. [45]. 

In this chapter, the computational aspects of the TNS zonal algorithm are 
presented, and then the mathematical and topological considerations of separated 
flows will be discussed to better understand them. After this, the impact of the 
artificial dissipation and turbulence models on separated flow will be demon- 
strated. Finally, the simulation of two low-aspect-ratio wings will be presented 
and comparisons will be made where data are present. 

4.2 Mathematical and Topological Considerations 

4.2.1 Background 

The character of a flow near a surface, to a large extent, can be inferred 
from the surface pattern or “skin-friction lines” which are the lines everywhere 
parallel to the wall shear stress vectors. The pattern of skin friction lines can 
often be approximately determined by the experimental oil-flow technique [85]. 
However, the thickness of the oil film plays an important role in the appearance 
of the resulting oil-flow pattern. If the film is thick and accumulates in certain 
regions of the surface, it may lead to misconceptions (private communication 
with Dean Chapman, 1984). 

Although the definition of separation and reattachment is well established 
in two-dimensional steady flow's by the existence of a reverse flow, a straight- 
forward extrapolation of two-dimensional concepts to three-dimensional flows is 
not possible. In three-dimensions, the meaning of the term “reverse” becomes 
ambigious. Our present understanding of three-dimensional separation stems 
essentially from flow-visualization experiments. If, for instance, a surface oil- 
flow technique is used, Tobak and Peake j 86] indicate that the convergence of 
oil-streaks onto a particular line is a necessary condition for three-dimensional 
separation. But, whether this is also a sufficient condition is a topic of current 
debate. 

Maskell [87 ] discussed the global structure of three-dimensional flows in qual- 
itative terms, calling the pattern of limiting streamlines (the streamlines very 
close to the surface which lie closely along the skin-friction lines) the “skeleton” 
of the entire flow field. A precise mathematical framework to conceptualize the 
pattern of “limiting streamlines” was proposed by Legendre [88] in which the 
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limiting streamlines were considered to possess the properties of a continuous 
vector field. One principal property is that only one limiting streamline can pass 
through a regular (that is, nonsingular) point in the flow field. On the basis 
of this hypothesis, the elementary singular points of the field can be categorized 
mathematically (see Sec. 4.2.2). Then, the types of singular points, their number 
and the rules governing the relations between them can be said to characterize 
the pattern. 


In this context, three-dimensional separation has been defined by the con- 
vergence of skin-friction lines onto a particular skin-friction line that originates 
from a singular point of particular type, the saddle point [85], [86]. However, 
this may not be the only possibility. Under given circumstances, it appears 
that a separation line may also emanate from a nodal-type singular point [86]. 
The issue was clarified more concisely by Lighthill [85j by tying the postulate 
of a continuous vector field to the pattern of skin-friction lines on the three- 
dimensional body rather than to the limiting streamlines just above the surface. 
There are two advantages to working with the skin-friction lines: 1) they are 
defined uniquely everywhere on the surface, even in the vicinity of lines of sep- 
aration, which are themselves skin-friction lines, and 2) with skin-friction lines 
being defined uniquely everywhere on the surface, the pattern of skin-friction 
lines can be viewed as a continuous vector field. 


4.2.2 Mathematical Foundations for the Skin-Friction Field 


In order to examine the conditions for separation and behaviour of critical 
points, it is necessary to consider the wall shear stress field. Skin-friction lines are 
defined as the integral curves of the wall shear stress vector exerted by the fluid on 
the wall. If at any point the wall shear stress vanishes, the field possesses a saddle 
or a nodal point, at which the skin-friction lines’ directions are indeterminate. 
The “phase-plane” and “phase-space” methods of exploring the properties of 
solutions of ordinary differential equations have proved to be extremely successful 
in the field of nonlinear dynamical systems. Oswatitsch [89], and Lighthill [85] 
examined viscous flow patterns close to a rigid boundary and classified certain 
classical points which can occur. They used a mathematics which is equivalent 
to the phase-plane trajectory analysis. Also Perry and Fairlie [90] and Hornung 
and Perry [91] applied these techniques to fluid flow problems. The following 
analysis outlines the mathematical procedure for the topography of skin-friction 
lines and critical points. In the analysis, cartesian coordinates are assumed for 
demonstration purposes, however the same conclusions could be drawn if the 
curvilinear coordinates were used. 
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The system of equations for particle paths in a steady flow is given by 


dx 

dt 


= u{x,y,z ) 


dy 

dt 

dz 

dt 


= v{x,y,z) 


(4.1) 


where x,y, and z are Cartesian coordinates and u,v, and w are the Cartesian 
velocity components along x, y, and z respectively (see Fig. 14). The skin-friction 
lines, or limiting streamlines, are obtained by exploring the limiting behavior of 
Eqs. (4.1). But because of the no-slip boundary condition at the wall, Eqs. (4.1) 
are trivial (i.e. u = v = w — 0, as z — » 0). At this point, it is helpful to introduce 


a pseduo-time dr such as 

II 

(4.2) 

Rewriting, Eqs. (4.1) yield 

dx u 



* dr z 

dy v 

y=T = - 

dr z 

dz w 

Z = d t = 'z 

(4.3) 


It is necessary to analyze the limiting behavior of the shear-stress tensor given in 
Eq. (3.3). At the surface, the following conditions for the velocity components 
and their derivatives hold: 


u — v — w — 0 

U z = Uy = V x = Vy = W x = Wy = 0 

but from the continuity equation, w z = 0 also. Hence Eqs. (3.3) reduce to 


(4.4) 
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du 

dz 

dv 

dz 


(4.5) 


Taking the limits as z — > 0 for the last two of the preceeding equations yields 
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or 


u 




(4.7| 


in which u — (u,v) T , t w = (r xz ,r yz ) T where the superscript T denotes the 
transpose of a matrix, and f{x,y) and g(x,y) are nonlinear functions. 

It is assumed that, for steady flow of an incompressible Newtonian fluid with 
constant viscosity at finite Reynolds number over a smooth continuous wall, the 
velocity and pressure are regular, so that local solutions can be written as Taylor 
series expansions in the space coordinates. The mathematical support for this 
assumption is discussed by Ladyzhenskaya [92]. Eq. (4.7) reveals that the vector 
quantity n/z is equal to T w /p, and the vector field u /z has the same integral 
curves as the wall shear stress; i.e., the skin-friction lines. When u / z is expanded 
into a Taylor series about the critical point, and only the lowest order terms are 
retained, the following linear system of first-order ordinary differential equations 
is obtained: 

r~n r„ u i r„i 

(4.8) 
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or 


F x 


(4.9) 

The matrix F of Eq. (4.9) has eigenvalues Aj and A 2 which may in general be 
real or complex. The corresponding eigenvector slopes are 


m x = (Ai - a) jb -- c /( Aj - d) 
m 2 - (A 2 - a)/b - c/(A 2 - d ) 


(4.10) 


which, for the case of Aj and A 2 real, may be shown to correspond to the slopes 
of certain trajectories which emanate from the critical points. 


There is just one skin-friction line through each point on a surface, except for 
a point of separation or attachment, where t w — 0. These points are the singular 
points of the differential equations governing the topography of the skin-friction 
lines. Such singular points are classified (see Kaplan [93]) depending on the 
values of the following quantities: 

p = — (a + d) = -trace F 
q (ad — be) = det F 

and 

Ai, 2 = ~\\PT (p 2 - 4 9 ) 1/2 ] 

The classification of possible critical points is shown in Fig. 15 from Ref. [90]. A 
singular point where q < 0 is a “saddle point”; a point where q > 0, however, is 


(4.11) 

(4.12) 
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a “nodal point.” A nodal point (Fig. 15a) is a point where all of the skin-friction 
lines except one (labeled AA in Fig. 15) are tangential to a single line BB. A 
focus (Fig. 15b) differs from a nodal point in that it has no common tangent 
line. An infinite number of lines spiral around the singular point, either away 
from it (a focus of attachment) or into it (a focus of separation). At a saddle 
point (Fig. 15c), only two particular lines, CC and DD, pass through the singular 
point. The flow directions on opposing sides of the singular point are inward on 
one particular line and outward on the other line. All of the other lines pass 
by the singular point in directions consistent with the directions of the adjacent 
particular lines. Critical points corresponding to values of p and q along the axes 
(p = 0 and/or q = 0) and on the parabola p = 4 q 2 are “degenerate” forms. 

4.2.3 Topography of Skin Friction Lines 

Singular points have certain characteristics that largely determine the dis- 
tribution of skin-friction lines on the surface. The nodal point of attachment is 
typically a stagnation point on the forward face of the body, where the free-stream 
attaches itself to the surface. Hence the nodal point of attachment behaves like 
a source from which the skin-friction lines issue and circumscribe the body. The 
nodal point of separation, however, behaves like a sink where the skin-friction 
lines on the body surface terminate. In other words, this is the rear stagnation 
point of the body. The saddle point acts typically to divide the skin-friction 
lines issuing from these nodes. The total number of singular points for a possible 
pattern on a smooth surface is subject to a topological rule that the number of 
nodal points must exceed the number of saddle points by 2 (see Kaplan [93]). A 
physical explanation for this is due to Lighthill [85]: since there are infinite num- 
ber of skin-friction lines on the surface, and they must begin and end somewhere, 
there is at least one nodal point of attachment and one nodal point of separation. 
If there are two nodal points of attachment, however, the skin-friction lines from 
each node must somewhere run into each other; this necessitates an introduction 
of a saddle point in between them. 

Figure 16 shows such a combination of two nodal points of attachment and 
one saddle point of separation. If there are n nodal points of attachment, there 
will be (n — 1) saddle points accordingly. Hence the whole combination is equiv- 
alent to a single node of attachment which behaves like a “source” as explained 
before. Similarly, m nodal points of separation and (m — 1) saddle points behave 
like a “sink” into which the skin-friction lines emanating from the “source” van- 
ish. Therefore, there are (m + n) nodal points and (ra + n — 2) saddle points, and 
the topological law is satisfied. A particular skin-friction line emerging from a 
saddle point prevents the lines from the nodal points of attachment from crossing 
each other. This particular skin-friction line is called a line of separation. Skin- 
friction lines from either side asymptotically converge on the line of separation. 


52 



Hence, it is possible to define the line of separation as a particular skin-friction 
line toward which other skin-friction lines converge rapidly (in an asymptotic 
manner, but not a cusp-like manner, c.f. (85], [90]); conversely, a line of reattach- 
ment is defined as a particular line from which other skin-friction lines diverge 
rapidly. It should be kept in mind that this definition is associated with the 
necessary condition of the separation. Whether it is also sufficient is a matter of 
current debate [86]. 

In the vicinity of the line of separation, the limiting streamline must leave 
the surface. This can be shown as follows: at a very small distance 2 from the 
surface the velocity was found in Eq. (4.7) 

u _ Tu, 

2 [l 


Continuity must be satisfied in a streamtube, for which the volume flow rate M 
is expressed as 

M = zd- — const. (413) 

2 

where d is the distance between the streamlines and q = | u |. From Eq. (4.7) 


where t w 


1 

Q- - ZT W 


t w |. From Eqs. (4.13) and (4.14) 
2 « (dr w ) _1/2 


(4.14) 


(4.15) 


Hence, the height of the limiting streamline 2 above the surface grows rapidly 
as the line of separation is approached. According to Eq. (4.15), the mechanism 
which causes the limiting streamlines to leave the surface is twofold: first, the 
resultant skin-friction t w becomes vanishingly small near the saddle point, and 
second, the distance d between the adjacent limiting streamlines, falls rapidly as 
the limiting streamlines converge toward the line of separation. 


In the remaining parts of this section, a series of attached and separated 
flow computations will be presented. The methodology of studying these flows 
will be from three points of view: 1) the actual phenomena (i.e., the experi- 
mental flow pictures), 2) the “conceptualization” or the “postulation” of these 
pictures in light of topological laws and rules, and 3) the computations. In any 
kind of interpretation, these three topics must support each other to provide a 
plausible explanation. The surface phenomenon is very important because the 
three-dimensional flow field above the surface depends on the surface flow. For 
instance, three-dimensional stream surfaces emerge from the separation lines, 
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whereas rolled-up vortex sheets emanate from spiral separation nodes, i.e. the 
foci on the surface. However, the three-dimensional separation phenomenon is 
very complex and a given surface pattern can produce many three-dimensional 
flow patterns depending upon interpretation. Dallmann [94] explains with exam- 
ples that there is no unique relationship between the pattern of wall streamlines 
and the flow field above the wall. He adds that topologically different flow fields 
may exhibit the same oil-flow picture. 

4.3 Effects of Artificial Dissipation and Turbulence Models 

Turbulence and artificial dissipation models play crucial roles in the simu- 
lation of separated flows and shock/boundary layer interactions. Since there is 
an intimate relation between the skin friction field and flow separation, correct 
computation of boundary layers is very important. In order to assess the relative 
effects of these models on the boundary layer characteristics, a series of test cases 
has been computed for attached and shock-induced separated flows. 

4.3.1 Effects of Artificial Dissipation 

The impacts of artificial dissipation schemes are firstly investigated for an 
attached airfoil flow. The test case consists of a two-dimensional flow around a 
NACA 0012 airfoil at = 0.5, Re - 2.89 x 10 c and a = 0°. The TNS program 
was run for a very large aspect ratio wing making the flow at the symmetry plane 
essentially two-dimensional. The computed and experimental pressure coefficent 
distributions are compared in Fig. 17. The experimental data is taken from 
Thibert et al. Ref. [95]. The agreement for this admittedly easy case is very 
good in terms of pressure. 

In Fig. 18, the turbulent boundary layer profiles computed by TNS, using 
different smoothing models and algorithms, are compared with those of an an- 
other CFD code, called TRIVIA. Ideally, the computed boundary layer profiles 
should be compared with the experimental ones. However, this is not an easy 
task because of the lack of experimental profile data. An alternative is to com- 
pare the computed profiles against those of another “proven” CFD code. For 
this purpose, the transonic airfoil interaction code of Steger and Van Dalsem 
[96], called “TRIVIA” was chosen. TRIVIA utilizes the same turbulence model 
as TNS and proved to be very accurate for a fairly broad range of airfoil flow 
conditions. 

The y + values of the first grid points off the surface were about 1 or 2 for 
TRIVIA and 2 or 4 for TNS. One remarkable feature of this comparison is that 
although all methods give the same C p distribution, the computed TNS boundary 
layer profiles using the diagonal algorithm of Section 3.3.2 are significantly less 
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full than either TRIVIA or TNS-block-algorithm boundary layer profiles. The 
TNS diagonal algorithm utilizes a nonlinear artificial dissipation model, while 
the block algorithm uses a constant coefficient artificial dissipation model. The 
discrepancy between the results may be explained as follows: the amount of 
artificial dissipation must be large enough to keep the algorithm stable, yet small 
enough not to excessively damp the flow physics. The nonlinear smoothing used 
in the diagonal algorithm was specifically designed for inviscid calculations [97]. 
Using this algorithm in Navier-Stokes calculations produces too much smoothing 
in the boundary layer. This is because the spectral radius o )y k,i used in Eq. 
(3.27) is in the order Af -1 , and since the normal grid spacing Af is very small 
in the viscous boundary layer, Oj t k,i becomes very large. This problem may be 
resolved by forcing the artificial viscosity to zero near solid boundaries. The 
nonlinear explicit smoothing terms of the diagonal algorithm given by Eq. (3.27) 
is modified to 


D.U = + (4.16) 


M) 


where 


f(M) 


! M 2 , for attached flows 
M, for separated flows 


(4.17) 


and M is the local Mach number. The suggested linear and quadratic variations 
of Mach number depending upon the type of flow was found experimentally. 
However, this reflects the fact that the amount of artificial dissipation in a sep- 
arated flow should be larger than its counterpart in an attached flow. This is to 
secure the numerical stability of the algorithm. Utilizing this concept proved to 
be very helpful, as evidenced in Fig. 19 where the modified dissipation model 
results compare favorably with the original TRIVIA results. 

Since the modified smoothing algorithm changes the boundary layer, it is 
expected that the skin-friction field, and consequently the separation character- 
istics, will also change. Therefore, another test case was set up to simulate a 
separated flow. The impact of the artificial dissipation, as well as the impact 
of the turbulence model, will be investigated on this flow. This case involves 
a transonic flow field about a NACA 0012 section wing with an aspect ratio of 
3.0, a taper ratio of 1.0, 20° of sweep, a Reynolds number based on chord of 8 
million, a free-stream Mach number of 0.826, and an angle of attack of 2°. The 
solid wind tunnel walls, which are modeled in this calculation (see Fig. 8), are 
slightly less than 2 chords above and below the wing. This case corresponds to 
the experimental work of Lockman and Seegmiller [61], and a comparison of the 
numerical results with this experiment follows. 
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The surface pressure coefficent distributions at three selected spanwise lo- 
cations are compared in Fig. 20. It is seen that the free-air shock position is in 
error by approximately 10% of chord, and inclusion of wind tunnel wall effects 
has shifted the shock downstream, in good agreement with the experiment. Also, 
the pressure plateau downstream of the shock, which signifies the shock induced 
separation, was not reproduced due to the turbulence model. Note also that 
general agreement between the wind tunnel simulation and experiment is better 
inboard of mid-scmispan than it is outboard. 

The emphasis of this work is on flow separation. In Fig. 21, an oil flow 
photograph taken from the experiment of Lockman and Seegmiller [61] shows 
the complexity of the phenomenon. In order to analyze and comprehend the 
separated flows, postulation of the skin friction lines is very helpful. The reader 
is reminded that Masked [87] called the skin-friction lines the “skeleton” of the 
entire flow field (see Sec. 4.2.1). The author’s depiction of the skin friction lines 
for the flow of Fig. 21 is presented in Fig. 22 to facilitate a possible correlation 
between the oil flow photograph and the numerical simulations. The oil-flow 
photograph displays a very interesting skin-friction map typical of transonic sep- 
arated flows. A strong swept shock wave extends almost over the entire span of 
the wing, and induces flow separation over two-thirds of the span. The swept 
shock wave is especially strong in the outer one-half of the span and causes a 
“mushroom-type” separation zone. In this zone, two counter-rotating vortices 
are separated by a saddle point, as necessitated by the rules for skin-friction line 
topography (see Sec. 4.2.3, and also Lighthill [85]). In addition to this interest- 
ing feature, two nodal points of separation, one nodal point of attachment and 
three other saddle points are recognized in Fig. 22. This skin-friction map seems 
to be plausible since no topological rule is violated. For instance, no saddle- 
point-to-saddle-point connection, which was shown to be structurally unstable 
[98] occurs, nor is there a nodal-point-to-nodal-point connection without a saddle 
point between them (see Fig. 16). 

The computed very-near-surface particle paths, or skin-friction lines, which 
resemble the oil flow pattern of the experiment are presented in Figs. 23 and 
24. The pattern in Fig. 23 was computed using the diagonal algorithm with 
the standard nonlinear smoothing. This simulation does compare well in terms 
of the extent of the separation region, which has been underpredicted. Indeed, 
the experimental separation is more than twice as large as the computed one. 
The pattern in Fig. 24, however, was produced using the same algorithm but 
with modified smoothing. The modified explicit nonlinear smoothing coefficient 
given by Eq. (4.16) was employed, and f(M) = M was used as suggested by Eq. 
(4.17). Note that the result of Fig. 24 is in better agreement with the experiment 
as indicated by the larger spanwise extent of the separation bubble. Although 
there is improvement, topologically Figs. 23 and 24 are the same; both have one 
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nodal point of separation and one saddle point. The nodal point is on the line 
of separation and the saddle point is on the line of reattachment. 

4.3.2 Effects of Turbulence Model 


The turbulence model also has important effects on the surface pattern. 
The Baldwin-Lomax model, which is an equilibrium model, predicts too sharp 
an increase in the outer layer eddy viscosity through the separation shock [99]. 
The increase of eddy viscosity from its local equilibrium value far from the wall 
leads to more turbulent mixing and to more shear stress to balance the adverse 
pressure gradient, suppressing the tendency toward separation. Rotta [100] con- 
cluded from experimental data that when turbulent flow is perturbed from its 
local equilibrium state, a distance of about one order of magnitude greater than 
the boundary layer thickness is required to attain a new equilibrium state. To 
account for the upstream turbulence history effects, Shang and Hankcy [101] 
used a relaxation eddy viscosity model, i.e., 


Mr = MT eq - (p-T eq - PT 0 ) ex P 



(4.18) 


which was applied on the y = const, lines, where is the turbulent dynamic 
eddy viscosity, Mr e «j ' s the local equilibrium eddy viscosity evaluated from Eq. 
(3.34), ht o * s the eddy viscosity at upstream location Xq, and A is the relaxation 
length. A good review of turbulence model relaxation techniques can be found in 
Hung [99]. Conceptually, Eq. (4.18) approximates the experimental observation 
that, in an abrupt disturbance of a turbulent flow, the Reynolds stress remains 
nearly frozen at its initial value while it is being convected along streamlines, and 
then exponentially approaches a new equilibrium state. In a numerical calcula- 
tion, the initial location of the disturbance from which the relaxation process is 
initiated, Xq, and a relaxation length scale which describes the exponential decay 
of the eddy viscosity distribution, A, must be specified. There are two limiting 
cases which bound the relaxation length. For A = 0, the turbulent eddy viscosity 
equals the local equilibrium value, and for A = oo, the initial value n t 0 is frozen 
and is used everywhere in the region x > x,,. 


The relaxation turbulence model has been applied to the separated flow case 
in the previous comparison. Varieties of A values were tried, and A = 40 6o was 
found to give a good agreement with the experimental pressure coefficient dis- 
tribution as well as with the skin-friction field. Here 6q is the thickness of the 
boundary layer at xq- The resulting computed skin-friction field for the case 
discussed in Figs. 21-24 is presented in Fig. 25. As seen from the figure, the 
streamwise separation is somewhat substantially increased; in fact, it is some- 
what overpredicted. Also the location of the separation bubble is more inboard 
than the experimental one, but this may be due to insufficient tip resolution in 
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the computational grid. The remarkable feature of Fig. 25, however, is that 
it reproduces the dominant feature of the experiment, that is, the mushroom- 
shaped separation zone. This feature was lacking in the Figs. 23 and 24. In 
Fig. 25, although the experiment shows two foci on each side of the saddle point, 
the computation predicts two nodal points of separation. This is topologically 
acceptable, although not accurate, because a node and a focus are topologically 
equivalent (see Fig. 15). The saddle point on the reattachment line was also 
predicted by the computation. 

The impact of the relaxation turbulence model on the pressure coefficient 
is displayed in Figs. 26 and 27. In Fig. 26, several sectional pressure coeffi- 
cients are shown across the separation zone. The solid line shows the relaxation 
turbulence model and dashed line shows the equilibrium turbulence model. Fig- 
ure 26 reveals that the equilibrium model predicts strong shock waves which 
more closely resemble the shocks captured in inviscid flow models. The relax- 
ation model, however, predicts the pressure plateau downstream of the shock 
due to the shock-induced separation. Another feature of the comparison is that 
the shock wave system of the relaxation model moved slightly forward relative 
to that of the equilibrium model. Nevertheless this is not totally unexpected. 
The displacement thickness <5* of the boundary layer increases due to the shock- 
induced separation. This portion of the boundary layer acts like a wedge which 
displaces the shock wave system forward. In Fig. 27, computations from the two 
turbulence models are compared with experiment. The computed shock posi- 
tions are somewhat forward of the experimental ones. However, it is remarkable 
that the pressure plateau predicted by the relaxation model is very close to the 
experimental one in Fig. 27b. On the other hand, the equilibrium model does 
not predict such a phenomenon. 

There are several things to keep in perspective when interpreting the dis- 
crepancy between the experiment and computation. One of the main problems 
is, of course, the turbulence modeling, as was proven by the improved prediction 
of the pressure plateau by using a simple relaxation model. Another significant 
factor in comparing the computation and experiment is that the experiment may 
need Mach number and angle of attack corrections, a common requirement when 
comparing with experimental data. Thirdly, the tip resolution in the TNS pro- 
gram is still poor which may be responsible for the inaccuracy of the location of 
the separation zone. Finally, it is known that prescribing wind-tunnel exit pres- 
sure in the computations is useful for predicting the extent of separation and the 
shock location [23]. Much improved results would be expected if these sources 
of error were eliminated. In the future, a more detailed study will be made with 
much finer grid resolution using more powerful machines together with a better 
turbulence model. Nevertheless, in this stage, it should be emphasized that the 
relaxation turbulence model was applied to the computation of three-dimensional 
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separated flows on wings for the first time, and improvements were observed in 
the skin-friction topography as well as the pressure coefficient of a low aspect 
ratio wing. 

In the remaining parts of this chapter, two low aspect ratio wings will be 

studied. Of these wings, the first one is an advanced technology wing called 
“WING C,” and the other wing is a rather simpler wing with NACA 0012 airfoil 
cross-sections. The flow topology on these wings will be thoroughly investigated 
by running many cases and comparing with the experimental data when they are 
available. The techniques developed in this section were not employed on these 
cases presently. However, in the future, as was explained in the previous para- 
graph, an attempt will be made to rerun these cases with higher grid resolution 
and better turbulence model. For the present, the aforementioned limitations 
should be kept in perspective when interpreting the results. 

4.4. Simulation of WING C Flow Fields 

4.4.1 WING C Design and Testing 

In recent years, there have been significant efforts to compute flow fields 
around wings in the transonic regime. Consequently, it is important to assess 
the accuracy of these computations by comparing them with reliable experimen- 
tal data. NASA Ames Research Center and the Lockheed-Georgia Company 
contributed to this effort by conducting a cooperative computational and exper- 
imental investigation of an advanced technology wing, called “WING C” (Fig. 
3). The designers of WING C, R. Hicks of Ames Research Center and B. L. 
Hinson of Lockheed-Georgia, used two existing computer codes: FL022 [102], 
a three-dimensional full potential equation solver and a numerical optimization 
program based on the Method of Feasible Directions [103]. A highly swept, low- 
aspect-ratio wing was selected that had supercritical airfoils with relatively thick 
sections, moderate aft loading, mild shock waves, and a mild pressure recovery. 
The mild shock wave was accomplished by limiting leading edge local Mach num- 
bers to maximum of 1.2 to the leading edge. The design conditions selected were 
a Mach number of 0.85 and a lift coefficient of about 0.5, at an angle of attack 
of 5°. 

A small-scale 0.261 meter semispan model of WING C was tested in the 
Lockheed-Georgia Wind Tunnel (CFWT), which has a (0.508 x 0.712m.) test 
section. The test Reynolds number based on the mean aerodynamic chord was 
10 million. Surface pressures were measured both on the wing and on wind- 
tunnel walls for comparison with calculations of wall effects. The small-scale 
model data are presented in Ref. [104], and a comparison of these measurements 
with several three-dimensional transonic inviscid codes is presented in Refs. [105] 
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and [106]. The top and side walls of the three-dimensional test section have a 
variable porosity capacity of up to 10%. Most of the wind tunnel testing was 
conducted at a fixed wall porosity of 4 percent for minimum wall interference 
effects. The influence of the wind tunnel walls on the test data was explored 
and a Mach number correction of A M = -0.005 was found necessary. An angle 
of attack correction of Act = 0.9° was also needed in the experiments to match 
the computations. The method of matching leading-edge pressures was used to 
select an experimental angle of attack of 5.9°, for which the experimental and 
computed (FL022 code, with a — 5°), leading-edge pressures agreed. Transition 
strips were located at a fixed distance from the leading edge equal to 5% of the 
mean aerodynamic chord on both the upper and lower surfaces of the wing. 

A series of tests with the WING C was also performed by NASA Ames 
Research Center as a part of the cooperative effort. For the Ames tests [107], 
a large-scale 0.90 meter semispan model was tested in a (1.8 x 1.8 m.) wind 
tunnel. Model blockage ratio in the test section was 1.3% at zero angle of attack. 
Surface pressure measurements and oil-flow studies were made at the design 
angle of attack of 5° over a Mach number range of 0.25 to 0.96 and a Reynolds 
number range of 3.4 x 10° to 10 x 10°. No Reynolds number effect on the results 
was reported to exist. The lift interference from the tunnel walls was reported 
to be small. This is because the leading-edge pressures of the experiment and 
computations in the correlations happened to agree with each other at the design 
angle of attack of 5°. Transition strips were installed at 4.5% chord. 

4.4.2 Attached Flow Cases 

The first set of WING C flow computations using the TNS program to be 
discussed will be attached flow cases. These cases are included to demonstrate 
the two-dimensional nature of the flow in the attached flow regime. This two- 
dimensional flow aspect was part of the design goal undertaken by the WING C 
design project. Attached flow calculations were done for two supercritical cases. 
In all these and subsequent computations, the Reynolds number is based on the 
mean aerodynamic chord of the WING C, and turbulent flow calculations are 
started at the very leading edge without any transition model. 

The first case consists of a slightly supercritical flow at M a 0 — 0.70, a = 5°, 
and a Reynolds number based on mean-aerodynamic-chord of 6.8 million. Figure 
28 shows the pressure coefficient distribution of WING C at 5 spanwise stations. 
The agreement between the two experiments and computation is quite good even 
at the r] — 0.90 (tip) location. However, the TNS results agree with the small- 
scale model data [104 ] slightly better than with the large-scale model data [107], 
The discrepancy between the two sets of data is probably caused by differences in 
wall interference effects, and differences associated with the measurement of angle 
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of attack and free-stream Mach number. An increase in the effective boundary- 
layer thickness decambers the aft section of a supercritical airfoil, which can 
significantly affect the predictions of lift and pitching moment. Figure 28 and 
the following figures showing pressure coefficient distributions reveal that the 
effect of the thickening boundary layer near the trailing-edge is more significant 
in the large-scale model data than it is in the small-scale model data. Therefore, 
the trailing-edge pressures of the large-scale model data are consistently lower 
than those of the small-scale model data. Although the general agreement is 
good, the TNS trailing-edge pressure values are overpredicted compared to the 
experiments. Figure 29 shows the experimental oil flow picture from Ref. [107] 
and the computed upper-surface skin-friction lines side by side. As seen from 
the figure, the flow is almost entirely two-dimensional on the wing, and the 
computation accurately reproduces this situation. 

The second computation is again an attached flow case with higher Mach 
number. The flow conditions are = 0.82, a - 5°, and Re m . a . c . = 6.8 x 10 6 . 
The computed pressure coefficient distributions are compared with data from 
the small-scale and large-scale model tests in Fig. 30. The comparison is very 
good in general, with the TNS solution in slightly better agreement with the 
small-scale model data [104]. The same considerations about the discrepancies 
between the experiments and computation presented for the M 0 0 = 0.70 case 
apply here. It is also useful to understand the shock wave development on the 
wing for interpreting shock-induced separations. Mach number contours, which 
show the shock wave pattern on the upper wing surface, are presented in Fig. 
31. A shock wave forms parallel to the leading-edge, and a small second shock 
wave forms nearly perpendicular to the root at the tip. Figure 32 shows an 
experimental oil-flow photograph and the computed skin-friction lines. In the 
oil-flow picture, only a weak shock wave is evident, indicated by a slight S-shape 
in the oil streak lines which lie in the outboard section of the wing between 
about 15% and 25% chord. This shock wave is not strong enough to separate 
the boundary layer, contrary to the earlier computations by Mansour [34], The 
flow is almost two-dimensional in the sense that the spanwise flow is negligible, 
and flow deflection angles are less than 10° except near the leading-edge. 

4.4.3 Separated Flow Cases 

The first separated flow case consists of the WING C design conditions: 
Moo = 0.85, a = 5°, and Re - 6.8 x 10 6 . These conditions were intended to 
result in attached flow with a mild shock wave and a mild pressure recovery. But, 
for reasons discussed in Ref. [107], these conditions produced a “local” (as called 
by its author) flow separation. In this dissertation, the adjectives such as “local” 
versus “global” or “open” versus “closed” will be avoided because of the current 
debate. However, the interested reader may refer to Ref. [86] for particulars. 
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The computed pressure coefficient distributions compared with experimental 
data are presented in Fig. 33. Two shock wave systems are distinguished: a 
leading edge shock and a secondary shock. The computed leading-edge shock 
agrees well with the experiments, though there are some discrepancies in the 
second shock. The computed second shock seems to agree best with the average 
of the experimental results in this region. Again the viscous effects decamber 
the aft portion of the wing and the trailing-edge pressures of the large-scale 
model data (Keener [107]) are lower than the small-scale model data (Hinson 
and Burdges [104]). Note that the TNS trailing-edge pressures agree well with 
the small-scale model data. Figure 34 shows the wing-planform Mach number 
contours plotted at a location above the upper-wing surface which displays the 
well-known transonic lambda shock-wave pattern. 

The experimental oil-flow picture of the same case is presented in Fig. 35. 
Figure 36 shows the postulated skin-friction field. A separation line caused by 
the swept shock wave emanates from a saddle point with two counter-rotating 
vortices on either side. The conjectured skin-friction map also features two saddle 
points and one nodal point of attachment. In three-dimensional separation, it is 
not possible to define a closed separation zone or “bubble,” as in two-dimensions. 
Here, the separation zone is largely fed by the vortical flow on the inboard side 
of the saddle point of separation. The tip flow is highly curved inboard towards 
the separation zone but does not enter the zone. Also, some streamlines are en- 
trapped by the inboard vortex while others pass by the separation region without 
being trapped. The flow is almost two-dimensional outside the separation zone. 

The computed skin-friction lines for this case are presented in Fig. 37. Note 
that an angle of attack correction of +0.9°, i.e. a = 5.9°, was used in the com- 
putations to get a closer computational/experimental agreement. In this com- 
putation, the global features of the experiment are predicted well: the location 
and size of the separation line, the streamlines being trapped by the vortex-like 
formation inboard of the separation line, the curvature of the tip streamlines and 
almost two-dimensional flow outside the separation zone. The critical points of 
the skin-friction map weren’t predicted, except a nodal point near the corner-like 
region between the line of separation and the tip flow streamlines. Overall, the 
simulation is globally very good and is quite encouraging. Moreover, although the 
critical points were not predicted, the computed topology is capable of producing 
a map similar to the experimental one: first of all, by increasing the grid resolu- 
tion and using the techniques described in Sec. 4.3, the vortical flow inboard of 
the line of separation can possibly evolve into a vortex and, since there is already 
a nodal point on the other side, emergence of a saddle point between them is 
inevitable. Such an evolvement was shown to be possible in Sec. 4.3 by applying 
new artificial dissipation techniques and turbulence models. These techniques 
weren’t applied in this case, but if applied, a significant chance for improvement 
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would be possible. Also, it should be kept in mind that the grid near the tip 
is coarse and increasing the resolution would probably improve the correlation. 
Another important issue to address is that the thin-layer approximation may 
be deficient in this region where there are strong streamwise and crosswise vis- 
cous flow gradients which are not accounted for in the thin-layer approximation. 
In addition, experimental uncertainties associated with wind-tunnel wall interfer- 
ence, angle of attack measurement, free-stream Mach number measurement, etc., 
must be considered in evaluating the computational/experimental correlations. 

In light of these encouragements, it is interesting to go one step further and 
study the flow field above the surface. Investigation of a three-dimensional flow is 
by nature very complex and needs sophisticated graphics utilities. In this work, 
a post-processing program called PLOT3D developed by Buning [108] at NASA- 
Ames was used. PLOT3D is a graphics program capable of plotting all relevant 
quantities of the flow field in three-dimensions. Particularly, the particle-tracing 
utility of the PLOT3D was helpful because of its power to display the actual 
fluid motion. 

The particle paths, in general, can be divided into two groups: “streamlines” 
and “sectional streamlines.” These two terms are explained in Fig. 38. If a par- 
ticle is released from point A, its path may be plotted in several ways, including 
the three shown, AB, AC, and AD. If its motion is not confined to any surface, 
it moves according to the fluid dynamic forces and becomes a “streamline.” Note 
that in steady flow, pathlines are equivalent to streamlines. If the motion of the 
particle is confined to x — y or y - z planes, “sectional streamlines” are obtained: 
the lines AD and AC respectively. It is not hard to imagine a sectional stream- 
line. For example, the line AD may be likened to the motion of a free-surface 
fluid particle when a body is immersed into a liquid and set into a motion. In 
this case, the motion of the particle is confined to the free-surface of the liquid. 
The sectional streamlines are computed from only the components of velocity on 
that particular surface, i.e., the projections of the velocity vectors on the sur- 
face. Note that the sectional streamlines are not the same as the projection of 
streamlines on the surfaces. For example, if the streamline AB was projected on 
the x — y or y — 2 surfaces, lines other than AC or AD would be obtained. 

Now, let us explore how to obtain various lines and surfaces which are of 
importance in the topography of flow fields. First of all, in steady flow, the 
trajectory of a “free” particle is a streamline. It shows exactly where the fluid 
is going. The skin-friction lines, which resemble the oil-flow, are simply the 
sectional streamlines traced on a surface very near a body. The reader is reminded 
that as z — ► 0 the u /z field becomes identical to the surface stress field of a body 
(see Sec. 4.2.2). Therefore, in Computational Fluid Dynamics, the skin-friction 
lines are obtained by plotting the particle paths constrained to the first grid 
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surface off a body. However, the distance 2 from the body must be essentially 
constant to make the u / z field identical to the u field. 

The physical meaning of the sectional streamlines is limited. The sectional 
streamlines coincide with the true streamlines only on the solid surfaces and 
symmetry surfaces if the flow is also symmetric. However, keeping in mind the 
limitations, the sectional streamlines are still instrumental for understanding the 
flow fields. At this point, some figures which will be discussed later in this paper 
can be used here as examples. As for example in Fig. 68, the sectional streamlines 
are used to illustrate the vertical extent of the separated region, and the vortical 
flow inside that region. From the very same figure, it could be imagined that 
a line which emerges from the line of separation extends above the wing, and 
divides the vortical flow in the separated region from the free stream which flies 
over that region. 

If we define a surface which emanates from a line of separation, and which di- 
vides the mass of separated fluid from the surrounding fluid as a “streamsurface,” 
the aforementioned line in Fig. 68 which is composed of closely packed stream- 
lines, is the cross-section of that “streamsurface” versus the JL plane. Therefore, 
we understand that it is possible to define the boundaries of a streamsurface by 
cutting it through with a sufficient number of planes each of which displays the 
sectional streamlines. This was done, for example, in Figs. 54 and 55. The sepa- 
rated zone could also be studied by using the sectional streamlines in a different 
orientation. In this case, as in Fig. 56, the sectional streamlines are plotted in 
the JK (body-parallel) planes. These figures serve to define approximately the 
boundary of the streamsurface as observed from above. 

In other examples, Figs. 62 and 67, an even more dramatic application of 
the sectional streamlines has been made. If we concentrate, say, on Fig. 67, the 
boundaries of the vortical flow are evident. In this figure, particles are released 
from different normal levels and are confined to stay on their respective planes. 
Had the Figs. 67 b, c, d been placed on top of each other, the boundaries of the 
vortical flow would have coincided. On the other hand, the true streamlines which 
were traced for the same case in Fig. 65 don’t give us such clear information, 
but give a very complex picture. Hence, the sectional streamlines are extremely 
useful for such cases. 

There are certain things to keep in mind when dealing with the sectional 
streamlines, however. Since they are the integrated curves of the projections 
of the velocity vectors on a certain surface, the orientation angle of the cutting 
plane is very important. Unless the cutting plane angle is nearly perpendicular 
to the lines of separation or reattachment, a distorted view is obtained. This is 
illustrated in Fig. 39. The figure depicts skin-friction lines for an ideal case: a 
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separation “bubble” which might occur on an infinitely swept wing. The skin- 
friction field at the root of such a “bubble” is cut by three planes I-I, II-II, and 
III-III, making an angle of 90°, -45°, and 45°, respectively, with the lines of 
separation and reattachment. The separation and reattachment points are des- 
ignated by S, S/ S, M and R,R/ R", respectively, on the three cutting planes. 
Observe that only plane I-I accurately represents the separation and reattach- 
ment lines (S-R). For the other planes, the plotting program cannot accurately 
predict these lines, showing the distance between the singular points to be shorter 
on plane II-II (as S'-R'), and longer on plane III-III (as S"-R"). The situation 
for the actual three-dimensional flows is even more difficult and, without use of 
other sectional streamlines, it is not possible to gather accurate information. 

These particle tracing techniques are very helpful in understanding three- 
dimensional flow fields, and they will now be employed to survey the separated 
WING C flow. Figure 40 shows the WING C surface grid together with the index 
designations for the planes of interest which will be used in subsequent figures. 
Figure 41 is the same as Fig. 37, except it is a three-dimensional perspective view. 
In this figure, only one particle at each spanwise r) station was released from the 
J = 15 grid line parallel to the leading-edge. The particles were confined to the first 
grid surface off the body, and therefore become sectional streamlines which define 
skin-friction lines. This figure defines the separation line and the particles which 
are entrapped by the separation zone, and also shows that the tip streamlines 
do not enter the zone. To help understand the following particle-path plots, 
the difference between the streamlines and sectional streamlines on the wing are 
demonstrated in Fig. 42. Instead of releasing just one particle from each spanwise 
station (as in Fig. 41), two particles were released in different modes. One set of 
particles was confined to the first grid surface off the body, and the second set 
was released without constraint from the same physical position. Hence, the first 
set becomes sectional streamlines or, to be more specific, skin-friction lines, and 
the second set becomes streamlines. As observed, two particles released at the 
same point trace the same path for a short distance and then diverge. If the two 
particles are released upstream of the separation zone, they trace similar paths 
until they encounter the separation. In the vicinity of separation, the constrained 
particle is entrapped but the free particle is not. It is easy to understand why 
the unconstrained particle flies over if the separation zone is thought of as a 
solid barrier. This situation may be likened to the phenomenon in Fig. 38. 
Here, line AD resembles a constrained particle path and line AB resembles an 
unconstrained particle path. Observe that the separation zone in Fig. 42 is thin, 
as deduced from the curvature of streamlines above the region. 

A series of pictures in Fig. 43 displays the separation surface emanating from 
the separation line as well as other interesting features of the three-dimensional 
boundary layer. In each picture, particles are released along an L-grid line which 
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is perpendicular to the wing from near the body to above the boundary layer. 
The particles are released from spanwise stations K = 14, 15, 18, and 21, which 
span the separation zone. It should be emphasized that these particles are con- 
strained to stay on their respective £-77 parallel surfaces to define the separation 
surface. If any particle hits the separation surface, it must go around it because 
the separation zone behaves like a solid object. Figure 43a shows the particles 
released from the K = 14 station, which is the last station, for which particles 
are not entrapped by the vortical-like flow of the separation zone. The invis- 
cid streamlines are almost parallel to the free-stream direction. As the viscous 
boundary layer is gradually penetrated, the streamline curvatures rapidly in- 
crease toward the separation zone. Figure 43b shows how the particles released 
from K = 15 go around the separation zone. The first particle path off the body 
is the skin-friction line which represents the line of separation. Notice that the 
second particle path follows parallel to the line of separation and becomes a tip 
flow. Figure 43c shows the particles released from the K = 18 plane. This fig- 
ure clearly shows how the separation surface evolves from the line of separation. 
When the particles released from different levels approach the separation* zone, 
they can not penetrate it; rather they flow over the separation surface, which is 
also a surface made of streamlines. Finally, Fig. 43d shows that the particles 
released from the K = 21 plane avoid the separation zone and pass by it without 
entrainment. This was already discovered in the study of the oil-flow pictures 
(see Figs. 35 and 36). To give an idea about the thickness of the separation 
zone, sectional streamlines constrained to lie within f f planes are presented in 
Fig. 44. The separation zone is extremely thin, being less than 0.001 chords. To 
complete the flow field survey, sectional streamlines in planes parallel to the wing 
surface (£-/?) are presented in Fig. 45. The separation and reattachment lines 
are observed in the L = 2 plane, which is the first grid surface off the body. The 
Figs. 45 b, c and d, define the boundary of the streamsurface w'hich emanates 
from the line of separation in Fig. 45 a. The separation line-type lines in Figs. 
45 b, c, and d are the cross-sections of the streamsurface with the JK planes. 


The second separated case associated with WING C consists of a flow at 
off-design conditions: — 0.90, a — 5°, and Re - 6.8 x 10 G . The computed 

and experimental pressure coefficient distributions are compared in Fig. 46. 
The leading-edge and aft shock-waves are stronger than those of the previous 
case, and the aft shock-wave is further downstream. This is a consequence of 
the increased free-stream Mach number. Despite a slight overprediction of the 
trailing-edge pressures, the agreement between the experiment and computation 
in this difficult case is very good except in the tip region, probably because of 
the coarse grid effect. The strength and location of the lambda shock-wave is 
accurately predicted. The planform view of WING C showing the lambda-shaped 
shock-wave is presented in Fig. 47. This agrees very well with the experimental 
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pattern reported in Ref. [104]. The Mach number contours of the (-<; plane 
which cuts the wing in the streamwise direction at the station 77 = 0 . 77 % is 
shown in Fig. 48. The zonal boundaries are also shown in this plot. Note 
the smoothness with which the shock wave crosses the zonal interface boundary. 
This indicates that the communication between the blocks is implemented in a 
conservative manner. In addition, most of the other contours cross the zonal 
interface boundaries in a smooth and continuous fashion. Downstream in the 
wake where the fine viscous grid interfaces with a relatively coarse grid (see Fig. 
5), the wake abruptly stops. This aspect of the solution exists because the coarse 
inviscid grid cannot retain the sharp wake gradient. 

The surface oil-flow photograph and the present postulation of the skin- 
friction lines for this case are presented in Fig. 49. When the Mach number 
is increased from 0.85 to 0.90, the flow separation that existed over the outer 
30% of the wing moves slightly downstream and grows significantly in size. The 
counter-rotating vortices which emerged in the previous case grow substantially, 
and the larger vortex inboard of the separation line extends as far as the trailing- 
edge. A line of separation emerges from the saddle point of separation. There 
is one focus on the tip side and a node of separation on the inboard side of this 
line. Another saddle point exists very near this node of separation. This saddle 
point lies between the node and the second, larger, focus. There exists another 
saddle point near the trailing-edge. The inboard vortex is strong and provides 
a mechanism to feed fluid to the separation zone by entrapping the streamlines 
coming from the inboard side of the separation zone. The flow between the 
separation zone and the symmetry plane is again almost two-dimensional in the 
sense that the spanwise flow is negligible. 

The computed skin-friction lines for the same case (A / 0 0 = 0.90) are pre- 
sented in Fig. 50. As observed, the agreement between these computed lines 
and the experimental picture in Fig. 49 is rather poor. A weak swirling motion 
is exhibited in the computation on the inboard side of the separation zone, yet 
the smaller vortex near the tip is nonexistent. Because of the coarseness of the 
present grid near the tip, this vortex is probably not resolvable. The disagree- 
ment occurs close to the tip where the computation does not capture the shock 
location well. It is this poor shock location (see Fig. 46) which moves the sepa- 
ration line too close to the trailing edge and causes the disagreement between the 
experiment and computation. In an attempt to better align the computed line 
of separation with that of the experiment, this case was recalculated at a slightly 
lower Mach number. It was found that when the Mach number was reduced 
slightly, to 0 . 88 , significant changes were observed which were in close agreement 
with the Moo — 0-90 experiment. 

The computed skin-friction lines for this slightly lower Mach number are 
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shown in Fig. 51. This flow pattern is indeed more descriptive of the experimental 
oil-flow at Mqo = 0.90: the line of separation moves forward, in better agreement 
with experiment, and there is substantially stronger vortical flow in the area of 
the inboard vortex. The node of separation in the experiment was captured. The 
saddle point which lies near the trailing-edge of the wing in the experiment was 
also captured. However, the vortex close to the tip has not been resolved. Had 
this vortex been resolved, emergence of the experimental saddle point between 
this vortex and the node would have been automatic according to the continuity 
equation. The computation also predicts a reattachment line at the trailing- 
edge near the tip, though this doesn’t exist in the experiment. To recapitulate, 
although quantitative details are not accurately reproduced many qualitative 
details are. The location and extent of separation are in reasonable agreement 
with experiment, as is the pressure distribution. The onset of vortex formation 
which is evident in the experiment is also present in the computation. In addition, 
experimental uncertainties associated with wind-tunnel wall interference, angle 
of attack measurement, free-stream Mach number measurement, etc., must be 
considered in evaluating the computational/experimental correlations. 


Although the surface flow was not accurately reproduced, it is instructive 
to explore the three-dimensional flow exhibited by the computation above the 
wing. Again, the sectional-streamlines will be used to study the flow field. The 
computed skin-friction lines which were presented in Fig. 51 are shown again in 
perspective view in Fig. 52 to give a more in-depth look at the flow field. The 
rapid convergence and divergence of the skin-friction lines clearly mark the sepa- 
ration and reattachment lines, respectively. To give an idea about the thickness 
of the separation, sectional streamlines are presented on a streamwise cut 
plane at spanwise station K = 17 (Fig. 53). The thickness of the separation zone 
is about 0.005 chords; which is much larger than the previous case. Figure 54 is 
an expanded view of the previous picture in the vicinity of the line of separation. 
Note that a shear surface lifts off the body from the line of separation at point 
S on the f-f computational plane. This shear surface partitions the outer flow 
from the separated flow, and is made of many streamsurfaces emerging near the 
separation line which enfold each other. In Fig. 55, a series of sectional stream- 
lines in (£-f) planes which span the separation zone are presented. As seen in 
the figure, the boundary layer displacement thickness grows from station K = 12 
to K - 20. In this figure, the vertical scale was expanded 5 times with respect to 
the horizontal scale to better observe the details. At the K = 19 and 20, a vortex 
core is observed under the shear surface. It should be kept in mind that, in all 
these plots, there is a substantial amount of transverse flow perpendicular to the 
f f planes. Also note that the points of separation and reattachment cannot be 
determined accurately from these views for the reasons enumerated in reference 
to Fig. 39. Accurate detection is only made possible by the skin-friction lines. 
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Sectional streamlines through the separation zone in the surface-parallel 
computational planes (f -?? surfaces) are presented in Fig. 56. This figure gives 
a sequence of views from the body surface to the edge of the separation zone. The 
shear surface emanating from the line of separation acts as a barrier between the 
flow within the zone and the flow upstream of it. In this figure, the approximate 
location of the shear surface is seen as a line of cross section between each body- 
parallel computational plane and the zone. 

The last computed case associated with WING C involved another off-design 
condition: Mqq = 0.95, with a — 5°, and Re = 6.8 x 10 6 . The comparison 
of the computed skin-friction lines with the experimental oil-flow photograph 
was similar to the M ^ = 0.90 case discussed before and is not presented here. 
However, the pressure coefficient comparison is shown in Fig. 57. Despite the 
high free-stream Mach number and the large amount of separation, reasonable 
agreement is obtained except at 90% semi-span. 

4.5 Simulation of NACA 0012 Wing Flow Cases 

Simulation of the NACA 0012 wing with wind tunnel walls modeled was pre- 
sented in Sec. 4.4. Experiments made by Lockman and Seegmiller [61] furnished 
the information necessary to assess the capabilities of the various features of the 
TNS code. Among these are the wind-tunnel wall modeling characteristics, and 
the effects of artificial viscosity and turbulence models. In order to ascertain the 
degree of robustness of the present algorithm and, in particular, the ability of 
the present zonal interface scheme to cope with large flow gradients, some NACA 
0012 wing free-air cases in more extreme conditions were run. 

Three of these cases are presented here: flows with free-stream Mach num- 
bers of 0.80, 0.85 and 0.90 at an angle of attack of 5 degrees with a Reynolds 
number based on chord of 80 million. Free-air conditions were necessary because 
the calculations with wind-tunnel walls modeled for the highest two Mach num- 
bers were choked (see Holst et al. [45]). As a result, no experiment is available 
to compare with these cases. However, very interesting separated surface flow 
patterns and three-dimensional flows above the wing were obtained. They are in- 
troduced here with the belief that these solutions will be helpful in understanding 
separated wing flows. 

The first case consists of a flow with mild streamwise separation. Flow 
conditions were M <*> = 0.80, o = 5°, and Re = 8 x 10 6 . Computed skin-friction 
lines are presented in Fig. 58. Except for significant inboard flow at the tip which 
is induced by the 5°angle of attack, the surface flow topology for this case is very 
similar to the wind-tunnel case presented in Fig. 23. Figure 58 describes a line 
of separation, a narrow streamwise separation zone, and a line of reattachment. 
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There is a node of separation and a saddle point on the lines of separation and 
reattachment, respectively. 

The second computation corresponds to a massively separated case. Flow 
conditions were = 0.85, ct = 5°, and Re — 8 x 10 6 . Figure 59 shows the 
wing upper-surface Mach number contours. As observed, there is an S-shaped 
swept shock wave which extends from the symmetry plane to very near the tip. 
This shock wave causes the boundary layer to separate between 20% to 90% of 
semi-span as shown in Fig. 60. A line of separation emanates from a saddle 
point with a focus on the inboard side and a nodal point of separation on the 
outboard side. This case exhibits a large reverse surface flow between the line of 
separation and the line of reattachment, which is located very close and parallel 
to the wing trailing-edge. Another saddle point can be recognized at the wing 
trailing-edge/tip. It seems that there is another line of secondary separation in 
addition to the primary streamwise separation. Surface flow coming from the 
tip region and from the line of reattachment coalesce on this line of secondary 
separation, and a streamsurface leaves the wing surface. 

The large vortex defined by the focus constitutes the root of a large vortical 
flow between the lines of separation and reattachment. Also, the streamsurface 
emanating from the secondary line of separation probably coils into two small 
secondary vortices above the wing. Development of this large vortex and the 
secondary vortices above the wing are illustrated in Figs. 61 and 62. Figure 61 
presents a series of three-dimensional pictures. They show sectional streamlines 
plotted on subsequent body-parallel computational planes in the JK plane. 

Pictures from a to h are successive planes above the body surface to a 
height of nearly 10% of a chord. The pictures a, 6, and c show lines similar 
to the separation line of Fig. 60. These lines are the cross sections of the 
computational planes with the separation surface which emerges from the line of 
separation. Note that there is a saddle point on each of these lines. The large 
vortex which coils up from the focus of Fig. 60 maintains its strength throughout 
the sequence of computational planes. The secondary separation line of Fig. 60, 
though not very clear, breaks into two secondary vortices as shown in Figs. 61d, 
61e, and 61 /. In these pictures, particles were released inside the separation zone 
to describe these vortices. Finally, pictures g and h of Fig. 61 show that these 
secondary vortices were probably dispersed and could not be traced any further, 
but the large vortex remains. Figure 62 complements Fig. 61 in describing the 
vortical flow pattern above the wing. In this figure, particles were released from 
the leading-edge at different L-levels and each was constrained in its respective 
JK plane. Figs. 62a, 62 b, and 62c show that the particles very near the wing 
are captured by the core of the large vortex, whereas the particles away from the 
body flow past it undisturbed. Figure 62 d shows the particles captured by the 
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big vortex and the secondary vortices as well. 

The last case is again a massively separated case. Flow conditions were 
Mqo = 0.90, a = 5°, and Re = 8 x 10 6 . Again in this case, an S-shaped 
shock wave extending from the symmetry plane to very near the tip causes the 
boundary layer to separate. This shock wave is stronger than the previous one, 
and separates the boundary layer in the same S-shape. In this case, computed 
skin-friction lines exhibit a very interesting surface pattern shown in Fig. 63. 
The separation zone extends from the symmetry plane to very near the tip, and 
the chordwise extent of the separation zone varies from about one- to three- 
quarters of a chord. A line of separation emanates from a saddle point which 
is located at about two-thirds of the span. A focus occurs near the tip, and a 
node of separation is located inboard at about 30% of semispan. On the line 
of reattachment, another saddle point is recognized. Expanded views of these 
critical points are shown in Fig. 64. In this figure, the focus quite interestingly 
defines a motion which is not swirling inward, but swirling outward from the 
core. Although not shown here, a surface pressure plot reveals that this vortex 
is located at a low pressure region. The primary mechanism which feeds this 
vortex seems to be the downward flow coming from the large vortex as seen in 
Fig. 65. In this figure, particles are injected from about 40% of chord at L = 
2 and are unrestricted. This figure exhibits two important phenomena: l) there 
are two vortices, a large and a small one, and they are located at either side of 
the saddle point of separation, 2) it seems as if the streamlines coming from the 
large vortex feed the small vortex. Note that particles by the small vortex eject 
from it at the tip side. This is a very curious mechanism. 

Unlike the particles in Fig. 65, the particles in Figs. 66a and 666 are 
restricted to their respective computational JK planes which are parallel to the 
surface. They are, therefore, sectional streamlines and define the core of the large 
vortex. The vortex core starts from the nodal point of separation shown in Fig. 
63 and expands as it moves out from the surface. Sectional streamlines are also 
used in Fig. 67 to give the situation a more three-dimensional appearance. In 
these figures, the same technique utilized in making Fig. 43 is used. Sectional 
streamlines on JK (body-parallel) computational planes are obtained by releasing 
particles from different levels at particular spanwise stations. In Figs. 67a and 
676, particles are released from the tip side of the saddle point at K = 18 and 14. 
Figure 67 a specifically shows that the particles closest to the surface are drawn 
into the small vortex. Particles more distant from the surface swirl around the 
small vortex but go to the core of the large vortex without entering the small 
vortex. Finally, particles away from the body run on straight paths without 
swirling. Figure 676 shows the same sort of behavior, but exhibits the particles 
captured by the large vortex more clearly. Figures 67 c and 67 d display the 
particle paths for particles released from K = 12 and 10, respectively, which are 
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inboard of the saddle point. 

The most important feature of these figures is that they define an approxi- 
mate envelope of the large vortex. The core of the vortex is observed at the bot- 
tom, and the near -circles at each level altogether define an approximate bound- 
ary of the three-dimensional vortex. This idea can be supported by placing Figs. 
676, 67c, and 67 d, on top of each other and observing that they fit together. Note 
that some particles are drawn toward the low pressure region of the small vortex. 
Finally, Fig. 68 shows the sectional streamlines in JL planes which are perpen- 
dicular to the wing surface. This cross-sectional plot was made at a semi-span 
station of 2y/b = 0.66. Note that at this station, the separation region is large 
and easily extends from the Navier-Stokes region, across the zonal boundary, 
into the Euler region. The particles pass smoothly across the interface boundary 
with no function or slope discontinuities. 




Fig. 14. The Cartesian coordinates, velocities and surface shear-stress vector 
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Fig. 15 . Classification of critical points. 
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Fig. 16 . A combination of two nodal points of attachment and a saddle point. 
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Fig. 17. Pressure coefficient comparison: NACA 0012 airfoil section, Ale — 0°, AR ^ 
TR = 1.0, Moo = 0.5, a = 0°, i?e c = 3 x 10 6 . 
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Fig. 19. Boundary layer velocity comparisons of TNS and TRIVIA: NACA 0012 airfoil 
section, Aj,# = 0°, large aspect ratio, TR = 1.0, = 0.5, a — 0°, Re c = 

3 x 10 6 . 
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Fig. 21. Oil-flow pattern on upper wing surface: NACA 0012 airfoil section, AR = 3, 
= 20°, &TW 1 ST 0°, TR = 1.0, A/qq = 0.828, a = 2°, Re c = 8 x 10 6 . 
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N s = NODE OF SEPARATION 
S = SADDLE 

N a = NODE OF ATTACHMENT 
F = FOCUS 


Fig. 22. “Postulated” skin friction lines for flow over the NACA 0012 wing of Fig. 21. 
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Fig. 25. Computed skin-friction lines for the flow over the NACA 0012 wing of Fig. 21 
with new dissipation model and relaxation turbulence model. 


COEFFICIENT OF PRESSURE 



Fig. 26. Comparison of computed pressure coefficents using the baseline and relax- 
ation turbulence models: NACA 0012 airfoil section, AR = 3, Ale = 20°, 
TR = 1.0, Moo = 0.826, a = 2°, Re c = 8 x 10 6 , a)2y/b = 0.51; b)2y/b = 
0.56; c) 2y/b = 0.62; d) 2yjb = 0.66. 
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Fig. 27. Comparison of experimental and computed pressure coefficients using differ- 
ent turbulence models: NACA 0012 airfoil section, AR = 3, A le = 20 , 
TR = 1.0, Moo = 0.826, a = 2°, Re c = 8 x 10 6 , a) 2y/b = 0.25; 6)2y/6 = 
0.50; c) 2y/6 = 0.78. 
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Fig. 28. Comparison of experimental and computed pressure coefficients for WING C: 
Moo = 0.70, (X — 5°, R^rn.a.c. ~ 0.8 X 10^. 


Fig. 29. An experimental oil-flow picture (from Ref. [107]) and computed upper sur 
face skin-friction lines for WING C: = 0.70, a — 2°, -Re m . a . c . — 6.8 x 10 6 
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Fig. 31. Mach number contours on the upper surface of WING C: Moo = 0.82, a — 
RCm.a.c. = 0.8 X 10°. 
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Fig. 32. Experimental oil-flow picture (from Ref. [107]) and computed upper surface 
skin-friction lines for WING C: = 0.82, a = 5°, Re m . a . c . = 6.8 x 10 6 . 
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Fig. 35. Experimental oil-flow picture (from Ref. (107)) and an expanded view of 
the counter-rotating vortices for WING C: Moo — 0.85, a = 5°, Re m .a.c. = 
6.8 x 10 6 . 
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Fig. 39. Illustration of obtaining different separation and reattachment points for a 
“bubble” depending upon the viewing angle. 
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Fig. 40 . Perspective view of the WING C surface grid and index designations for the 
computational planes. 
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. 42. Perspective view of computed “streamlines” and “sectional streamlines” for 
WING C: Moo = 0.85, a = 5.9°, Re m , a , c _ = 6.8 x 10 6 . 
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Fig. 43. 

i 


Sectional streamlines on 
“releasing” particles at 


JK (body-parallel) computational planes obtained by 
different L levels for a fixed value of J (jV/qq :=: 0.85). 
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Fig. 44. Sectional streamlines on JL (cross-sectional) computational planes obtained 
by “releasing” particles at different L levels from streamwise stations through 
the separation zone (M^ = 0.85). (Note the expanded z-scale.) 
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(b) L = 3 
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Fig. 45. Sectional streamlines on JK (body-parallel) computational planes. Particles 
are “released” from a line parallel to the leading-edge (Moo ~ 0.85). 







Fig. 46. Comparison of experimental and computed pressure coefficients for WING C 
Moo = 0.90, a — 5°, Re m . a . c . = 6.8 x I0 6 . 










Fig. 48. Cross-sectional Mach number contours for WING C at 2 y/b = 0.77: Af^, 
0.90, a = 5°, Re m . a . c . = 6.8 x 10 6 . 
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Fig. 53. Sectional streamlines at the K = 17 (cross-sectional) computational plane ob- 
tained by “releasing” particles at different levels from a station upstream of 
separation (M^ — 0.88). 
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Fig. 54. Expanded view of Fig. 53 near the separation point (z-axis has been expanded 
10 times more than the x-axis for better definition of the separation region). 
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Fig. 56. Sectional streamlines on JK (body-parallel) computational planes. The par- 
ticles are “released” from a line parallel to the leading-edge (M 0 0 = 0.88). 
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Fig. 56. concluded. 
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Fig. 60. Computed skin-friction lines for the flow over the NACA 0012 wing; M< 
= 5°, Re c = 8 x 10 G . 
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Fig. 61. Sectional streamlines on JK (body-parallel) computational planes. Particles 
are “released” from a iine parallel to the leading-edge (M'oo = 0.85). 
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Fig. 61. concluded. 
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Fig. 62. Sectional streamlines on JK (body-parallel) computational planes obtained by 
“releasing” particles at different L levels for a fixed value of J (Moo = 0.85). 
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Fig. 66. Sectional streamlines on JK (body-parallel) computational planes. Parti- 
cles are “released” from a line parallel to the leading-edge at 40% chord 
(Moo = 0.85). 
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Fig. 67. Sectional streamlines on JK (body-parallel) computational planes obtained by 
“releasing” particles at different L levels for a fixed value of J ( M ^ = 0.90). 
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CHAPTER 5 


CONCLUSION 


5.1 Summary 

The significant achievements in this work are associated with both numer- 
ical methods and fluid mechanics. A fast, efficient, and reliable computer pro- 
gram (TNS) was developed, and using this program, numerical solutions of the 
Euler/Navier-Stokes equations for the transonic flow about wing geometries were 
obtained. For the first time, computations produced good agreement with exper- 
iment for low-aspect-ratio, separated wing flows in the transonic regime. Under 
experimental guidance, new numerical techniques were developed to improve the 
computations. Contributions to the understanding of separated flows were made 
through the extensive use of computer graphics. Familiar as well as new fluid 
flow phenomena, which should be of interest to experimental and theoretical fluid 
mechanics, was displayed. 

The TNS computer code represents the first three-dimensional Euler/Navier- 
Stokes zonal algorithm. Using this program, fine grid-resolution flow fields about 
wing geometries have been achieved. It is a generic program which includes wind- 
tunnel wall or free air simulation capability, and is suitable, with specific exten- 
sions, for computing wing-fuselage or wing-fuselage-tail configurations. The TNS 
program is more than one order of magnitude faster than comparable Navier- 
Stokes solvers for wing geometries reported in the literature. It was shown that, 
instead of using the original block-ADI algorithm, by using a diagonal-ADI al- 
gorithm with certain enhancements, a speed-up of 40 was possible. 

Also, an efficient data management scheme using SSD (solid state device) 
resulted in an I/O time of about 1-2% of the CPU time for a typical run. On 
the other hand, if standard rotating disk had been used, an I/O time of about 
1000% of the CPU time would have resulted! 

Using the TNS program, good correlation between experiment and compu- 
tation was achieved for complex flow calculations associated with shock-induced 
separated flows. First, the effects of different artificial dissipation models incorpo- 
rated into the block and diagonal algorithms associated with the TNS code were 
assessed for attached and separated flows. It was found that these algorithms 
computed correct pressure distributions, but surprisingly predicted different skin- 
friction fields. Accurate computation of pressure is not the ultimate goal of a 
Navier-Stokes solver; even a good potential flow solver can give a reasonable an- 
swer if the flow conditions are not extreme. A Navier-Stokes code should resolve 
the boundary layer and reasonably compute the viscous and heat transfer effects. 
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In this dissertation, a simple scheme modifying the artificial dissipation model 
associated with the diagonal-ADI algorithm in the TNS code was proposed and 
successfully implemented for both attached and separated flow cases. As a result 
of this, for example, C j was increased roughly 50% for an attached flow case, 
and showed excellent agreement with the result of another “proven” viscous- 
inviscid interaction program. Also, using this scheme, improvements have been 
observed on the skin-friction maps of shock-induced separated flows as compared 
to oil-flow photographs. 

Another aspect of this work has been to prove that a simple idea, such as 
the relaxation turbulence model, can improve numerical computations. In this 
respect, certain limitations of the algebraic turbulence model due to Baldwin 
and Lomax associated with the TNS code were eliminated. As a result, the quite 
challenging task of obtaining pressure plateaus associated with boundary layer 
separation was achieved. Also, and perhaps more importantly, the distribution 
of critical points associated with the skin-friction map of one transonic separated 
case (NACA 0012 wing, M^ = 0.826, a = 2°) was changed from an incorrect 
topology into the correct topology as supported by the experimental oil-flow 
picture. 

However, the most important outcome of this work was the beginning of 
the fulfillment of a long-awaited expectation of seeing agreement between exper- 
imental oil-flow pictures and computations for separated flows. This desire was 
expressed by George Schairer, a retired senior Boeing company engineer, to the 
first author’s research advisor Brian Cantwell in a private communication. In 
this work, one very good correlation between the computed skin-friction field 
and experimental oil-flow for a highly swept, tapered, and cambered wing, i.e. 
the WING C, was obtained for the first time (in the case of Moo = 0.85, a = 5°). 
Also, two other computed cases (WING C, Moo = 0.88, a = 5° and NACA 0012, 
Mqo = 0.826, a = 2°) show relatively good correlation. All computed cases are 
in good agreement with experiment in terms of the pressure distributions. 

Investigation of the TNS code’s ability to compute massively separated flow 
cases in the transonic regime for the NACA 0012 wing resulted in some quite 
fascinating phenomena. The skin-friction maps of these cases produced very 
well-formed critical points, particulary two surface vortices; one spiraling inward, 
another spiraling outward! This is the first time computations of this type have 
been produced. It was learned (by private communication, Tobak and Fairlie, 
1985) that the presence of an outward spiraling vortex has been presumed to 
be less likely by quite a few experimentalists in their work. Also, it is believed 
that the illustrations of three-dimensional flow fields above the NACA 0012 wing 
and WING C via computed “streamlines” and “sectional-streamlines” put some 
light on the nature of three-dimensional vortical flows emanating from separated 
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turbulent boundary layers. This is an area where the experiments are limited to 
only surface oil-flow techniques, or again limited dye or smoke flow techniques. 
It is hoped that a theoretical or experimental fluid mechanics researcher may 
further benefit from the computed results. 

5.2 Recommendations for Future Work 

Throughout this work, the lack of sufficient grid resolution has been the 
single most important limiting item. Without enough grid resolution, it is not 
possible to give a definitive answer to the question of discrepancy between an 
experiment and computation. Therefore, the first recommendation is to rerun 
these cases with higher grid resolution when a faster computer with more memory 
is available. Another way of refining the grid without waiting for a new computer 
is to place additional zones in critical regions of the flow; for example, the tip. 
An H-C or O-C grid topology at the tip would be ideal for this task. However, 
this necessitates a more sophisticated interpolation program to handle the new 
situation. 

As might be remembered, another important shortcoming was the poor res- 
olution of the viscous wake. The viscous grid zones do not extend downstream 
far enough to resolve the far wake. Since the coarse Euler grid cannot support 
the high gradients associated with a wake structure, this causes certain inac- 
curacies, and can become a very significant factor especially at high angles of 
attack involving future applications with multiple lifting surfaces. However, it is 
expected that placing another zone in the wake will solve the problem. 

The importance of artificial dissipation was previously emphasized, and a 
simple correction was suggested which worked quite nicely. However, a better 
scheme is highly desirable. Probably, the present scheme which uses the sum of 
the spectral radii of the flux jacobians as a coefficient is not good, because it 
varies like 0(1/ Af) with distance f through the boundary layer. This scheme 
was originally designed for the Euler equations where it works well. An improved 
scheme which is especially designed for Navier-Stokes computations would be 
appropriate. 

Turbulence modeling is, of course, extremely important for high Reynolds 
number flows and for separated flows. It is known that, with the exception of 
few cases, the Baldwin-Lomax or Cebeci-Smith algebraic eddy viscosity models 
cannot properly handle two- or three-dimensional separated flows. Work needs to 
be done in this area implementing existing or developing new turbulence models. 
This should be done with sufficiently refined grids so as to eliminate coarse grid 
effects from affecting the results. 
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One subject which deserves special attention is wind-tunnel wall simula- 
tion. It would be appropriate to model porous wind-tunnel walls and assess their 
effects. Before this can be successful, detailed measurements at the walls are 
required. There is plenty of room for improvement with solid wind-tunnel walls 
as well. As is known quite well from two-dimensional computational work, pre- 
scribing wind-tunnel exit pressure as well as upper and lower wall pressures as 
boundary conditions is extremely important. It is believed that providing these 
capabilities to the TNS code should improve wind tunnel wall simulations. It 
should especially lessen the discrepancy between the experimental measurement 
and computation of the shock strength and location. 

All of the separated flow computations in this dissertation belong to shock- 
induced separations at high Reynolds and Mach numbers. However, it would be 
interesting to investigate flows for the subcritical or even low subsonic regimes 
at high angles of attack. Such a study would help in understanding separation 
due to adverse pressure gradients. The aircraft industry would benefit greatly 
from this, because such a calculation simulates stall behavior of a wing during 
take-off or landing. 
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Appendix A 


FLUX JACOBIAN AND VISCOUS COEFFICIENT MATRICES 

The flux jacobians A,B , and C are given by 


K o K x K 2 

K X (j > 2 — uO Kq + 9 — K 1(7 — 2 )u K 2 u — (7 — \)K\v 
K 2 (j > 2 - v 0 K x v - A 2 (7 - l)u K 0 + 6 - K 2 (i - 2 )v 

K 3 <t > 2 - wO K x w — A 3 (t i — l)u A 2 u> - K 3 { 7 - l)v 

9[2<^ - -,(e/p)J - * 2 | {if 2 [i(e/,,) - 

-(7-l)ud) -(7-l)i>0} 


K 3 0 

K 3 u - (7 — l)K x w K 1(7-1) 
A3V — (7 — \) K 2 w K 2 { 7 — 1) 

Kq + 0 - K 3 ( 7 — 2)w ^3(7 — 1) 

{A 3 [7(e/p) - </> 2 ] K 0 + 7# 
— (7 - l)u;0} 

where 

< f> 2 = 0.5(7 - 1)(« 2 + V 2 + w 2 ) 

9 = Xju + K 2 v + K 3 w 
and, for example to obtain A, 


•^0 — Aj — f x , K 2 = £ y , A 3 = 


also, the viscous coefficient matrix is given by 

' 0 0 0 0 0 - 

m 2 i ai$ f (p -1 ) a2^(p _1 ) a 3 < 5 ? (p -1 ) 0 

M n = m 3 1 a 2 <5 ? (/)' 1 ) as^jp -1 ) 0 (A.2) 

m 4 1 a 3 6^(p~ 1 ) a 5 <5 ? (p _1 ) a 6 ^(p _1 ) 0 

.m5i m5 2 m 53 ra 54 cto^Jp -1 ). 
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with 


m 2 1 = ai 6 { [-u/p) + a 2 6^{-v/p) + a z 6^{-w/p) 
m 3 i = a 2 «5 ? (-u/p) + a 4 6 ( (-v/p) + a 5 6^{-w/p) 
m 4 1 = a z 6f(—u/p) + a 5 6^(-v/p) + a 6 6^(-w/p) 
m 5 1 = ai6s(-u 2 /p) + a 2 6^(-2uv/p) + a 3 6^(-2uw/p) 

+ a 4 S^(-v 2 /p) + a 6 6^(-w 2 / p) + a 5 6 c (-2vw / p) 

+ oto6((—e/ p 2 ) + o:o^[(« 2 + v 2 + u> 2 )/p] 

(A3) 

m52 = — m 2 1 - ao 6 ^( u / p ), 17153 = —”^31 — a o ^<;( v / p ) 

m 54 = -m 41 - a 0 S^(w/p) 
a 0 = 7/cPr _1 (f 2 + £ 2 + d)>«i = m[(4/3)c 2 + 

<*2 = (A*/3) CasCy » a 3 = (m/3)?x£z 

= m[?* + (4/3)c 2 + ? 2 ], a 5 = (m/3 )?»?* 

a 6 = M[?x + Cy + (4/3)? 2 ] 
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Appendix B 


SIMILARITY TRANSFORMATION MATRICES 


The similarity transformation matrix for the Jacobian matrices A, B, and 


C is 


T k = 


h x 
k x u 

k x v + k z p 
k x w — k y p 

M>7(' 7-1) 


k y u - k z p 


kyV 


kyW + k x p 

\ky(f > 2 / (7 - 1) 


k Z U + kyP 

k z v — k x p 

k z w 

\kAVh - 1) 


+p(k z v — A: y w)] +p(k x w - k z u )] +p(k y u — k x v )] 


a 

a(u + k x c) 
a(v + k y c) 
a{w + k z c ) 


a 


a(u — k x c) 
a(t> — k y c) 
a(w — k z c ) 


(B. 1 ) 


a\(4> 2 + c 2 )/ (7 - 1) + c0\ a[((j> 2 + c 2 )/(7 - 1) - c0 ] J 

where, k = £, 77, and f for A, B, and C respectively. Also, the inverse transfor- 


mation matrix is given by 



’ [k x {l - <f> 2 /c 2 ) 
-p~ l {k z v - A v «;)] 

1M 1.- 0 2 A 2 ) 

-p~ 1 {k x w - k z u)\ 

(M 1 - <i> 2 lf) 

-p~ 1 (k y u - k x v)} 
0(<l> 2 - c0) 

. P(4> 2 + cO) 


k x {l - l)«c 2 
\~'k z p~ l 

+ ky{n - 1)«C“ 2 ] 

_ \-kyP~ 1 
+k z ( 7 - l)uc- 2 ] 
0\ k x c - (7 - l)u] 
-0\k x c + (7 + l)u] 


\k z p~ l 

+Mt ~ l)t>c~ 2 ] 
k y ( 7 - l)rc~ 2 

[k x P~ l 

+ k z {'i - l)vc -2 ] 
( 3 [kyC - (7 - l)v] 
-/ 3 [kyC+ (7 - l)t>] 


\~kyp~ 1 + k x ( 7 - l)wc~ 2 ] 
\kxp~ 1 + k y ( 7 - l)«;c _2 ] 
k z ( 7 — l)w;c -2 
- (7- l)w] 
-f3[k z c + (7 - l)u>] 


-Mt “ l)c -2 ' 

~ky{l ~ 1)C“ 2 

~k z { 7 - l)c -2 

- 1) 

- 1) 


(B. 2 ) 
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where, d — k x u + k y v + k z w and, for example, k x = k x /(kl + k 2 + k \ ) 1 / 2 , etc. 
Also, 4> 2 — 0 . 5(7 ~ l)(u 2 + v 2 + w 2 ), a = pjy/2c and /? = 1 /y/2pc. 
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